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Abstract 

Mark-recapture  analysis  of  populations  is  becoming  an  important  tool  in  population  biology. 
Mark-recapture  methods  can  be  used  to  estimate  transition  probabilities  among  life-stages 
from  capture  histories  of  marked  individuals  for  which  stages  can  be  determined  at  each 
sampling  occasion.  This  method  is  called  a  multi-stage  mark-recapture  (MSMR)  method. 
In  this  thesis,  I  describe  advances  I  made  in  the  MSMR  method  and  present  analyses  that 
apply  this  method  to  actual  data. 

The  advances  I  made  in  the  MSMR  method  are  motivated  by  a  need  to  provide  a  link 
between  mark-recapture  data  and  demographic  models  such  as  matrix  population  models 
and  integrodifference  models.  I  resolve  some  issues  that  are  commonly  encountered  during 
sampling,  such  as  the  fact  that  the  sex  or  life-stage  of  some  individuals  is  unknown  during 
some  sampling  occasions  and  that  individuals  become  unobservable  during  some  life-stages. 
I  introduce  a  stage-structure  that  permits  simple  conversion  of  estimated  transition  proba¬ 
bilities  into  a  matrix  population  model.  I  describe  an  algorithm  to  simplify  programming 
for  parameter  estimation.  I  also  introduce  a  method  to  estimate  the  distribution  of  dispersal 
displacements  (a  dispersal  kernel)  from  mark-recapture  data. 

I  apply  some  of  the  methods  described  above  to  data  of  the  North  Atlantic  right  whale 
(Eubalaena  glacialis ).  The  right  whales  are  considered  one  of  the  most  endangered  mam¬ 
mals.  The  current  population  size  is  about  300  in  the  northwestern  Atlantic,  and  the 
number  is  declining.  I  applied  the  multi-stage  mark-recapture  statistics  to  the  17-year  in¬ 
dividual  sighting  history  data.  Using  the  estimated  transition  probabilities,  I  constructed 
a  population  projection  matrix,  which  was  used  for  further  demographic  analyses.  I  found 
that  the  population  was  slowly  increasing  in  1980,  but  it  started  to  decline  slowly  around 
1992.  I  show  that  (1)  this  change  was  caused  by  increased  mortality  of  females  that  have 
just  given  birth,  (2)  protecting  two  females  a  year  from  the  deaths  is  enough  to  prevent  the 
declining  trend,  and  (3)  demographic  stochasticity  is  a  more  important  factor  influencing 
their  long-term  viability  than  environmental  stochasticity. 
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Chapter  1 

Introduction 


Mark-recapture  (capture-recapture)  methods  have  become  an  important  tool  in  population 
biology.  The  methods  are  used  to  estimate  population  size  (e.g.  Petersen,  1896;  Lincoln, 
1930;  Jackson,  1933),  survival  probability  (e.g.  Burnham  et  al.,  1987;  Lebreton  et  al.,  1992; 
Caswell  et  al.,  1999),  and  transition  probabilities  among  spatial  locations  (e.g.  Arnason, 
1972,  1973)  and  life  stages  (e.g.  Nichols  et  al.,  1992b).  In  a  mark- recapture  study,  a  popu¬ 
lation  consisting  of  previously  marked  (or  identified)  individuals  is  resampled  in  subsequent 
sampling  occasions.  By  resampling  the  population,  information  on  the  capturability  of 
individuals  is  obtained.  This  information  allows  better  estimates  of  population  parame¬ 
ters  than  methods  where  the  parameters  are  calculated  directly  from  count  data  under  the 
assumption  that  probability  of  capturing  individuals  is  1. 


1.1  Brief  History  of  Mark-Recapture  Methodology 

Mark-recapture  sampling  was  originally  developed  as  an  alternative  approach  to  a  census 
for  estimating  population  size.  In  the  original  method,  the  resampling  was  only  done  once 
to  obtain  the  numbers  of  individuals  that  were  marked  and  unmarked  at  the  resampling  oc¬ 
casion,  from  which  capture  probability  was  estimated.  By  dividing  the  number  of  originally 
marked  individuals  by  the  capture  probability,  the  abundance  of  individuals  in  a  population 
was  estimated.  Let  x,  y,  and  z  be  the  number  of  individuals  marked  at  time  1,  the  number 
of  individuals  caught  at  time  2,  and  the  number  of  recaptured  marked  individuals  at  time 
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2,  respectively.  Then  the  population  size  N  can  be  estimated  as 


N  =  —x,  (1.1) 

z 

where  |  is  an  estimated  capture  probability  assuming  it  was  the  same  at  the  two  sampling 
occasions.  Equation  (1.1)  is  often  called  the  Lincoln-Petersen  estimator  in  mark-recapture 
literature.  This  estimator  implicitly  assumes  that  the  population  is  closed  to  immigration, 
emigration,  birth,  and  death,  and  also  assumes  that  capturability  of  marked  and  unmarked 
individuals  are  the  same.  This  method  was  first  used  by  Laplace  in  the  18th  century 
to  estimate  the  human  population  in  France  (Cormack,  1968).  It  was  rediscovered  later 
by  Petersen  (1896)  and  Lincoln  (1930)  and  used  to  estimate  population  size  of  fish  and 
waterfowl,  respectively. 

The  method  to  estimate  survival  probability  from  mark-recapture  data  was  developed  by 
Cormack  (1964),  Jolly  (1965),  and  Seber  (1965).  This  method  is  often  called  the  Cormack- 
Jolly-Seber  (CJS)  method.  While  the  only  data  required  for  the  Lincoln-Petersen  estimator 
are  number  of  individuals  marked,  number  of  marked  individuals  recaptured,  and  the  to¬ 
tal  number  of  individuals  recaptured,  the  CJS  method  requires  additional  information  on 
whether  or  not  each  uniquely  marked  individual  was  recaptured  at  each  sampling  occasion 
(i.e.  individual  capture  history  data).  In  the  CJS  method,  the  individual  capture  histories 
are  modeled  assuming  all  marked  individuals  have  the  same  capture  probability;  however, 
it  does  not  make  any  assumption  on  the  capturability  of  unmarked  individuals.  In  reality, 
the  capture  probability  may  be  different  among  individuals  (i.e.  heterogeneous),  and  ones 
with  lower  capture  probability  may  be  more  likely  to  be  unmarked  than  others.  Because 
the  CJS  method  does  not  make  any  assumption  on  the  capturability  of  unmarked  individu¬ 
als,  it  is  more  robust  to  the  heterogeneity  in  capture  probability  than  the  Lincoln-Petersen 
estimator.  Extensive  reviews  of  the  CJS  method  can  be  found  in  Burnham  et  al.  (1987) 
and  Lebreton  et  al.  (1992). 

As  an  extension  of  the  CJS  method,  Arnason  (1972)  and  Arnason  (1973)  developed 
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a  method  to  estimate  movement  rates  of  individuals  among  discrete  spaces,  along  with 
their  survival  and  capture  probabilities.  This  method  is  called  multi-type  mark-recapture 
method.  The  multi- type  mark- recapture  method  requires  information  on  where  individuals 
were  marked  and  recaptured,  in  addition  to  the  individual  capture  histories  that  are  required 
for  the  CJS  method.  The  data  are  modeled  with  capture  probability,  survival  probability, 
and  transition  probabilities  among  spatial  locations.  This  method  was  later  generalized  by 
Nichols  et  al.  (1992b)  to  be  used  with  multi-stage  mark-recapture  data,  in  which  “stage” 
can  be  defined  in  terms  of  the  life  cycle  (age  classes,  size  classes,  developmental  stages),  or 
any  other  property  that  divides  individuals  into  subgroups. 

1.2  Overview  of  Thesis 

The  focus  of  my  thesis  work  is  on  advancing  the  mark-recapture  analysis  method  to  extract 
useful  information  for  demographic  analyses  from  the  data.  The  thesis  consists  of  five  main 
chapters  following  this  introduction.  The  first  three  chapters  are  primarily  on  new  mark- 
recapture  methodology,  and  the  last  two  chapters  are  application  of  the  methods  to  actual 
data  of  the  North  Atlantic  right  whale. 

In  Chapter  2,  I  present  four  major  improvements  in  the  multi-stage  mark-recapture 
method.  (1)  I  use  a  Markov  chain  formulation  of  the  life  cycle  to  express  the  likelihood 
functions  in  matrix  form,  which  makes  numerical  calculations  simpler.  (2)  I  introduce  a 
method  to  incorporate  capture  histories  with  uncertain  stage  and  sex  identifications,  which 
allows  the  use  of  capture  history  data  with  incomplete  information.  (3)  I  introduce  a 
simple  function  that  allows  multinomial  transition  probabilities  to  be  written  as  functions 
of  covariates  (time  or  environmental  factors).  (4)  Finally,  I  show  how  to  convert  transition 
probabilities  estimated  by  the  MSMR  method  into  a  matrix  population  model. 

In  Chapter  3,  I  present  a  new  method  to  estimate  survival  probability  from  the  capture 
histories  of  marked  individuals  that  become  unobservable  during  some  life  stages.  During 
mark-recapture  studies  of  open  populations,  animals  often  temporarily  emigrate  from  study 
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areas.  Such  temporary  emigration  can  cause  biased  estimates  of  survival  probabilities  and 
other  parameters.  The  new  method  uses  stage-structured  models  that  include  one  or  more 
stages  representing  the  individuals  that  have  temporarily  emigrated.  Although  not  all  pa¬ 
rameters  can  be  estimated  in  such  stage  structures,  some  important  parameters  are  still 
estimable.  The  estimability  of  parameter  values  is  determined  from  the  rank  of  the  Jaco¬ 
bian  of  the  likelihood  function.  I  apply  the  temporary-emigration  mark- recapture  method 
to  artificial  data,  representing  various  life  histories,  and  demonstrate  consistency  between 
actual  and  estimated  values.  The  method  presented  in  this  chapter  will  be  especially  useful 
for  studies  of  seabird,  sea  turtle,  and  marine  mammal  populations  where  breeding  is  delayed 
and  where  individuals  are  sampled  only  on  their  breeding  grounds. 

In  Chapter  4,  I  present  a  new  method  for  estimating  a  distribution  of  dispersal  dis¬ 
placements  (a  dispersal  kernel)  from  mark-recapture  data.  One  conventional  method  of 
calculating  the  dispersal  kernel  assumes  that  displacements  are  normally  distributed  (e.g. 
resulting  from  a  diffusion  process),  and  that  individuals  remain  within  sampled  areas.  The 
first  assumption  prohibits  an  analysis  of  dispersal  data  that  are  not  normally  distributed  (a 
common  situation);  the  second  assumption  leads  to  underestimation  of  dispersal  distance 
because  individuals  that  disperse  outside  of  sampling  areas  are  never  recaptured.  The  new 
method  eliminates  these  two  assumptions.  This  method  uses  integrodifference  equations  to 
express  the  probability  of  spatial  mark-recapture  data;  associated  dispersal  and  recapture 
parameters  are  then  estimated  using  a  maximum  likelihood  method.  The  accuracy  of  model 
selection  was  examined  using  Akaike  Information  Criteria  ( AIC) .  The  method  presented  in 
this  chapter  will  allow  re-examination  of  existing  data  sets  and  suggests  designs  for  future 
mark-recapture  experiments. 

In  Chapter  5,  I  apply  the  method  presented  in  Chapter  2  to  the  data  of  the  North 
Atlantic  right  whale.  The  North  Atlantic  right  whale  ( Eubalaena  glacialis )  is  considered 
the  most  endangered  large  whale  species.  The  current  population  size  is  about  300  in 
the  northwestern  Atlantic,  and  the  number  is  declining.  I  applied  the  multi-stage  mark- 
recapture  statistics  to  the  17-year  individual  sighting  history  data  collected  by  New  England 
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Aquarium,  MA.  Using  the  estimated  transition  probabilities,  I  constructed  a  population 
projection  matrix.  In  this  chapter,  I  also  present  various  demographic  analyses  based  on 
the  projection  matrix. 

In  Chapter  6,  I  extend  the  analysis  in  Chapter  5.  In  particular,  I  analyze  the  impor¬ 
tance  of  demographic  and  environmental  stochasticity  to  the  population  viability.  I  apply 
the  branching  process  calculations  (e.g.  Caswell,  2001),  calculations  based  on  a  diffusion  ap¬ 
proximation  to  the  changes  in  log  of  population  size  (e.g.  Tuljapurkar,  1990;  Dennis  et  ah, 
1991),  and  simulation  (e.g.  Caswell,  2001)  to  calculate  the  expected  time  to  extinction  from 
the  estimated  parameters. 

At  the  time  of  writing  this  thesis,  Chapter  2  and  3  have  been  accepted  to  Ecology. 
Chapter  4  is  to  be  submitted  to  Biometrics.  Chapter  5  has  been  published  in  Nature. 
Chapter  6  is  to  be  submitted  to  a  journal  that  is  to  be  determined. 
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Chapter  2 


Estimating  population  projection 
matrices  from  multi-stage 
mark-recapture  data 

2.1  Introduction 

Mark-recapture  estimates  of  survival  probability  have  been  applied  to  many  animal  popula¬ 
tions  (e.g.,  Lebreton  et  al„  1992;  Forsman  et  al.,  1996;  Weimerskirch  et  al.,  1997;  Hastings 
&;  Testa,  1998;  Caswell  et  al.,  1999;  Pease  &  Mattson,  1999),  and  this  method  has  become 
an  important  tool  in  population  management.  Mark- recapture  estimates  are  based  on  cap¬ 
ture  histories  of  individually  identified  animals,  which  contain  information  on  whether  or 
not  each  individual  was  captured  at  each  sampling  occasion.  For  example,  capture  history 
data  may  be  obtained  by  annual  observations  of  banded  birds  or  photographically  identified 
whales.  When  such  data  are  available,  mark-recapture  statistics  are  considered  one  of  the 
best  approaches  for  estimation  of  survival  probability. 

Modern  demographic  analysis  goes  beyond  calculating  survival,  by  breaking  the  life 
cycle  into  stages  (which  may  be  based  on  age,  size,  developmental  or  behavioral  states, 
physiological  condition,  spatial  location,  or  any  other  property  that  divides  individuals  into 
subgroups).  The  fate  of  individuals  is  described  in  terms  of  transition  probabilities  among 
these  stages,  and  those  transition  probabilities  form  the  basis  for  matrix  population  models 
(Caswell,  2001).  Nichols  et  al.  (1992b)  introduced  a  method  to  estimate  transition  probabil- 


ities  among  stages  from  mark-recapture  data,  which  we  call  the  multi-stage  mark-recapture 
(MSMR)  method.  This  method  extends  the  method  originally  developed  to  estimate  prob¬ 
abilities  of  movement  among  spatial  locations  (Arnason,  1972,  1973;  Brownie  et  ah,  1993; 
Lebreton,  1995).  For  the  MSMR  method,  in  addition  to  information  on  whether  or  not  each 
individual  was  captured,  the  capture  history  data  must  also  include  the  stage  of  captured 
individuals  at  each  capture  occasion.  MSMR  models  account  for  inter- group  heterogeneity 
in  survival  and  capture  probability  by  grouping  similar  individuals  into  stages.  The  devel¬ 
opment  from  single-stage  to  multi-stage  mark- recapture  statistics  parallels  the  development 
from  unstructured  to  structured  population  models.  In  fact,  one  motivation  for  the  statisti¬ 
cal  development  was  the  need  to  estimate  parameters  in  stage-structured  matrix  population 
models  from  mark- recapture  data  (Nichols  et  al.,  1992b). 

The  analysis  is  based  on  maximization  of  a  likelihood  function  that  depends  on  all  the 
possible  sequences  of  stage  transitions  compatible  with  an  observed  capture  history.  There 
can  be  very  many  of  these  sequences,  and  one  of  the  most  complicated  parts  of  the  method 
of  Nichols  et  ah  (1992b)  is  writing  them  all  down  with  their  associated  probabilities.  In  this 
paper,  we  describe  the  life  cycle  as  a  Markov  chain,  and  take  advantage  of  this  description 
to  write  the  likelihood  in  a  simple  matrix  notation.  A  sketch  of  this  method  was  given 
in  Caswell  (2001,  Section  6. 1.2. 2).  Here  we  give  a  complete  presentation,  and  extend  the 
method  to  incorporate  uncertainty  in  stage  and  sex  identifications,  which  allows  the  use  of 
capture  histories  containing  incomplete  information.  We  also  introduce  a  simple  function 
that  allows  multinomial  transition  probabilities  to  be  written  as  a  functions  of  covariates 
(e.g.,  environmental  variables  or  time).  Finally,  we  show  how  to  convert  the  estimated 
transition  probabilities  into  a  matrix  population  model. 

2.2  MSMR  Statistics 

The  MSMR  method  involves  three  main  steps:  (1)  constructing  an  appropriate  stage  struc¬ 
ture,  (2)  expressing  the  likelihood  function  in  terms  of  parameters,  based  on  available  cap- 
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ture  histories,  and  (3)  finding  the  best  parameter  estimates  using  maximum  likelihood 
theory.  The  parameters  in  the  MSMR  model  are  those  that  define  the  capture  probabili¬ 
ties  of  each  stage  at  each  sampling  occasion  and  the  transition  probabilities  among  stages 
between  consecutive  sampling  occasions.  The  method  assumes  that  individuals  in  the  same 
stage  are  identical  and  independent,  but  that  individuals  in  different  stages  may  differ  in 
their  transition  and  capture  probabilities.  The  MSMR  method  is  very  flexible  and  can 
be  applied  to  almost  any  stage  structure.  Constructing  a  useful  stage  structure  that  is 
compatible  with  the  life  cycle  of  populations  requires  experience  in  addition  to  sufficient 
mathematical  and  biological  knowledge,  and  different  stage  structures  are  extensively  re¬ 
viewed  in  Caswell  (2001).  In  this  paper,  methods  for  expressing  the  likelihood  function 
and  estimating  parameters  are  described  assuming  an  appropriate  stage  structure  has  been 
constructed. 

To  make  our  discussion  more  concrete,  we  will  demonstrate  the  method  using  a  stage 
structure  (Fig.  2.1)  developed  to  describe  the  life  history  of  the  North  Atlantic  right  whale 
( Eubalaena  glacialis ).  This  is  a  two-sex,  multi-stage  model  that  distinguishes  calves  (stage 
1),  immature  individuals  (stage  2),  and  mature  individuals  (stage  3).  In  addition  to  these 
three  stages,  females  also  have  a  stage  for  individuals  nursing  a  calf  (stage  4);  we  call  the 
individuals  in  this  stage  “mothers.”  Stage  0  corresponds  to  death,  and  the  probabilities  (f>oi 
associated  with  the  arrows  going  to  stage  0  are  stage-specific  mortality  rates.  As  usual, 
“mortality”  includes  both  death  and  permanent  emigration. 

The  objective  of  the  MSMR  approach  is  to  estimate  the  transition  probabilities  associ¬ 
ated  with  each  arrow  and  the  capture  probabilities  of  each  stage.  In  the  next  section,  we 
show  how  to  construct  matrices  containing  the  transition  and  capture  probabilities,  and  to 
account  for  uncertainty  in  the  assignment  of  individuals  to  stages.  Then  we  will  show  how 
to  calculate  the  likelihood  in  terms  of  these  matrices. 


Female 


(b) 


Figure  2.1:  A  stage  structure  for  female  (a)  and  male  (b)  right  whales.  This  structure  is 
used  as  an  example  for  the  MSMR  statistics. 
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2.2.1  Transition  and  Capture  Probability  Matrices 

The  transition  matrix  is  constructed  by  first  putting  the  transition  probability  <f>ji  from 
(living)  stage  i  to  (living)  stage  j  in  the  (j,i)  position.  To  this  matrix  is  appended  a  row 
containing  the  probabilities  of  transition  from  each  stage  to  stage  0  (death)  and  a  column 
containing  the  probabilities  of  transition  from  stage  0  to  each  stage.  Because  we  treat  death 
as  a  stage,  the  result  is  the  transition  matrix  of  an  absorbing  Markov  chain,  with  death  as 
an  absorbing  state.  The  matrix  is  column  stochastic.  The  ability  to  treat  transitions  as  a 
Markov  chain  is  critical  to  our  analysis.  The  transition  matrix  for  females,  corresponding 
to  the  stage  structure  in  Fig.  2.1  is 


where  4>ji{t)  is  the  probability  of  females  making  transition  from  stage  i  to  j  between  time 
t  and  t  +  1.  The  upper  left  block  of  the  matrix  describes  the  transition  among  live  stages, 
and  lower  left  block  of  the  matrix  contains  stage  specific  probabilities  of  death.  The  1  in 
the  (5,5)  entry  is  the  probability  of  dead  individuals  remaining  dead  in  the  following  year. 
The  notation  <f>ji(t)  in  this  paper  corresponds  to  (f>tij  in  Nichols  et  al.  (1992b). 

Similarly,  the  transition  matrix  for  males,  corresponding  to  the  stage  structure  in 
Fig.  2.1,  is 
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We  have  written  and  as  a  function  of  a  vector  of  parameters  6.  Then  para¬ 
meters  could  be  the  (f>ji  themselves,  or  lower  level  parameters  from  which  the  <f>ji  can  be 
calculated.  The  objective  is  to  estimate  9. 

Capture  probability  matrices  Pt,  defined  for  females  and  males  separately,  contain  stage- 
specific  capture  probabilities  on  the  diagonal  and  zeros  elsewhere: 

(  Pi  (t)  0  0  0  0  N 

0  P2(t)  o  0  0 

p  [f\e)  =00  p3(t)  0  0  (2.3) 

0  0  0  p4(t)  0 

I  0  0  0  0  0; 

pi(f)  0  0  0 

0  P2(i)  0  0 

0  0  P3{t)  0 

0  0  0  0 

where  pj(t )  is  the  probability  of  capturing  individuals  in  stage  j  at  time  t.  This  notation 
corresponds  to  pjt  in  Nichols  et  al.  (1992b).  As  written  here,  P^  and  p|,”')  assume  that 
dead  individuals  are  never  captured  (p^  =  p q7^  =  0),  but  such  captures  could  be  included. 

Transition  and  capture  probability  matrices  can  be  defined  separately  for  females  and 
males  when  information  on  sex  identification  is  available,  as  in  our  example,  but  it  is  not 
always  possible  or  necessary  to  have  a  two-sex  model.  In  such  cases,  only  a  single  transition 
and  capture  probability  matrix  is  needed. 

2.2.2  Stage- Assignment  Matrices 

A  stage- assignment  matrix  is  defined  for  each  individual  each  time  it  is  captured.  The 
diagonal  elements  of  the  matrix  are  proportional  to  the  certainty  of  stage  identification  at 
time  t  (i.e.,  to  the  probability  that  the  individual  is  in  a  given  stage  when  it  is  captured). 
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This  probability  should  be  known  prior  to  estimating  transition  and  capture  probabilities. 
In  our  example,  individual  A;  is  a  female;  its  stage-assignment  matrix  is 

(  i4fc)(t)  0  0  0  o  N 

0  u(k\t)  0  0  0 

u{fe)  =00  u(k\t)  0  0  (2.5) 

0  0  0  u{k\t)  0 

,0  0  0  0  4fc)(t)  , 

'  » 

where  u^k\t)  is  the  probability  that  individual  k  at  time  t  is  in  stage  j  ( j  =  1,2, 3, 4). 
Similarly,  if  individual  k  is  male,  its  stage-assignment  matrix  is 

u{k\t)  0  0  0 

0  4fe)(i)  0  0 

0  0  u{k\t)  0 

ooo  4k)(t ) 

Because  we  assume  that  the  capture  probability  of  dead  individuals  is  zero,  the  value  for 
Uqk\t)  will  not  enter  into  the  likelihood  calculations.  Multiplication  of  Ujfe  by  a  scalar  has 
no  effect  on  the  maximum  likelihood  estimates. 

If  the  stage  of  the  individual  is  known  with  certainty,  its  stage-assignment  matrix  con¬ 
tains  a  1  in  the  corresponding  diagonal  entry  and  zeros  elsewhere.  On  the  other  hand,  if 

ffc) 

the  stage  of  an  individual  is  completely  unknown,  the  identity  matrix  can  be  used  for  U)  . 
This  specifies  a  uniform  probability  distribution  over  the  possible  stages.  Alternatively,  if 
an  independent  assessment  of  the  probability  is  available,  it  can  be  entered  into  the  matrix. 
For  example,  in  an  age-structured  model  of  fish,  age  of  fish  is  sometimes  determined  from 
their  length  using  age-length  keys  (e.g.  Fournier  Sz  Archibald,  1982;  Deriso  et  al.,  1985; 
Quinn  &  Deriso,  1999).  Such  a  key  could  provide  the  probability  distribution  of  age  of  the 
fish,  which  can  be  entered  into  the  stage  assignment  matrices,  for  an  age-structured  model. 
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2.2.3  Likelihood 


The  likelihood  of  the  parameter  vector  0  contains  contributions  from  the  capture  history  of 
each  individual.  We  denote  by  lk{0)  the  contribution  to  the  likelihood  from  individual  k] 
it  is  proportional  to  the  probability  of  the  capture  history.  That  probability  is  the  sum  of 
the  probabilities  of  all  possible  sequences  of  transitions  that  could  have  been  taken  by  the 
individual  k.  There  may  be  many  such  possibilities.  Their  sum,  however,  can  be  calculated 
using  the  transition,  capture  probability,  and  stage-assignment  matrices  by  the  following 
algorithm.  We  assume  the  individual  is  first  captured  at  time  t\. 

1.  Categorize  individual  k  by  its  stage  at  its  first  capture,  taking  uncertainty  in  stage 
assignment  into  account, 

U^e.  (2.7) 

where  e  is  a  vector  of  ones.  This  product  is  a  vector  whose  entries  are  proportional 
to  the  probabilities  of  the  initial  stage  of  the  individual  at  t. 

2.  Calculate  the  probability  distribution  of  the  stage  at  t2  by  multiplying  this  vector  by 
the  transition  matrix  . 

ttiUfe.  (2.8) 

3.  Calculate  the  probabilities  of  observation  outcomes  at  f2: 

if  individual  k  was  captured  at  t2,  multiply  by  the  sighting  matrix  P t2- 

PfaSt.U^e,  (2.9) 

if  individual  k  was  not  captured  at  f2,  multiply  by  (I  —  P t2): 

(I-Pt2)$flUt(f}e,  (2.10) 


where  I  is  the  identity  matrix, 
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4.  Account  for  stage  identification  at  t2  by  multiplying  by  the  stage  assignment  matrix; 
if  individual  k  was  captured  at  t2 

Ujfp^ufe,  (2.11) 

if  individual  k  was  not  captured  at  t2 

I(l-Pt2)*tiu£}e.  (2.12) 

5.  Repeat  steps  2-4  until  the  end  of  the  capture  history  for  individual  k.  The  result  is  a 
vector  whose  ith  entry  is  proportional  to  the  probability  of  all  the  pathways  by  which 
individual  k  could  have  moved  from  its  initial  stage  at  t\  to  stage  i  at  th  and  that  are 
compatible  with  its  capture  history. 

6.  The  final  step  is  to  sum  the  resulting  vector  of  probabilities  to  obtain, 

**(*)*  =  eTug}  (2.13) 

In  this  algorithm,  the  probability  distribution  of  the  individual’s  stage  is  updated  sequen¬ 
tially  over  time  taking  into  account  the  new  data  available  at  each  time  step  and  possible 
stage  transitions  determined  by  the  stage  structure.  Therefore,  the  right-hand  side  of  equa¬ 
tion  (2.13)  is  the  probability  of  the  capture  history  for  individual  k,  taking  into  account  all 
possible  transition  sequences  compatible  with  that  history. 

The  likelihood  lk(&)  is  calculated  using  only  female  or  male  matrices  if  the  sex  of  indi¬ 
vidual  k  is  known.  If  the  sex  of  individual  k  is  uncertain,  the  above  algorithm  is  repeated  to 
get  likelihoods  and  1^(0)  using  the  female-  and  male-specific  matrices,  respectively. 

Then,  the  likelihood  lk(&)  is 

lk(0)=pfl\P(9)  +  (1  -Pf)l{r\0),  (2-14) 
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where  pf  is  the  probability  that  the  individual  k  is  female.  The  probability  p/  is  1  or  0 
when  the  sex  of  the  individual  is  known  to  be  female  or  male,  respectively.  If  the  sex  of  the 
individual  is  unknown,  a  probability  must  be  provided  to  calculate  the  likelihood. 

Some  examples  of  probabilities  of  the  capture  histories  of  individuals  with  four  capture 
periods  are  shown  in  Table  2.1.  Because  our  example  contains  multiple  stages,  many  possible 
capture  histories  exist,  of  which  only  a  few  are  shown  in  Table  2.1.  For  simplicity,  we 
assume  that  sex  of  all  individuals  is  known  to  be  female.  None  of  the  likelihoods  in  Table 
2.1  contains  Pi,  because  the  probability  of  a  capture  history  is  always  conditional  on  the 
first  capture;  therefore,  capture  probability  at  the  first  sampling  time  cannot  be  estimated. 
For  the  same  reason,  the  likelihoods  of  individuals  5,  7,  8,  11,  12,  and  13  do  not  begin  at 
time  t  =  1,  because  capture  histories  prior  to  the  first  capture  of  an  individual  do  not  enter 
into  probability  calculations. 

Given  the  likelihood  functions  lk(9 )  for  all  individuals,  the  likelihood  associated  with 
the  data  consisting  of  n  capture  histories  is  proportional  to  the  product  of  the  n  likelihood 
functions 

L(0)  oc  J]  lk{6).  (2.15) 

k=  1 

Here,  we  assume  that  individuals  are  captured  and  make  stage  transitions  independently, 
but  based  on  identical  probability  distributions  (i.e.  we  assume  the  number  of  outcomes 
falling  into  the  possible  capture  history  sequences  is  multinomial). 

Maximum  likelihood  estimates  0  are  found  by  maximizing  L(6).  The  likelihood  function 
can  be  maximized  numerically  using  software  such  as  MATLAB.  For  example,  the  MATLAB 
routine  fminu  can  be  used  to  find  the  maximum  likelihood  by  minimizing  -  log  L[9). 

2.3  Transition  Probabilities  as  Functions  of  Covariates 

Transition  probabilities  <pji{t)  may  change  over  the  course  of  a  study,  and  the  changes  may 
be  correlated  with  various  factors.  We  would  like  to  model  the  probabilities  as  functions 
of  covariates  measuring  those  factors.  For  example,  population  density  and  sampling  effort 
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Table  2.1:  Some  possible  capture  histories  corresponding  to  the  example  stage  structure  in 
Fig.  2.1  and  their  likelihood. 


ID  (fc) 

Capture  History* 

Likelihood  lk{0) 

1 

3433 

eTui1)P4^3U^)P3^2U^P2^iU(1lje 

2 

122X 

eT(I  -  P4)^3U|2)P3^2U(2)P2#1ui2)e 

3 

22X3 

eTui3)P4^3(I  -  P3)*2u£s)P2*iU[3)e 

4 

1X22 

eTui4)P4^3u(4)P3^2(I  -  P2)*iUS4)e 

5 

X343 

eTui5)P4^3U$5)P3^2U(5)e 

6 

12XX 

eT(I  -  P4)#3(I  -  P3)^2U(6)P2^iUS6)e 

7 

XX33 

eTuf)P4#3U^7)e 

8 

X43X 

eT(I  -  P4)*3U?) P3$2uf  e 

9 

4XX3 

eTuf  P4^3(I  -  P3)*2(I  -  P2)^iUS9)e 

10 

2X2X 

eT(I  -  P4)*3l410)P3*2(I  -  P2)^iUS10)e 

11 

X3X4 

eTuin)P4#3(I  -  P3)^2U^n)e 

12 

XXIX 

eT(I-P4)#3U<12)e 

13 

X3XX 

eT(I-P4)$3(I-P3)^2U(13)e 

14 

2XXX 

eT(I  -  P4)*3(I  -  P3)*2(I  -  P2)^iUS14)e 

- - - 1 ; - rn~. - ~ 

Note:  When  stage  of  captured  individual  is  i,  Uj  is  a  matrix 

with  1  in  the  ith  row  of  the  ith  column  and  0  elsewhere. 


*  X  indicates  the  individual  was  not  captured;  numbers  indicate 
the  stage  of  captured  individuals. 
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were  used  to  model  the  survival  and  capture  probabilities  in  studies  of  the  roe  deer  ( Capre- 
olus  capreolus )  and  the  common  lizard  ( Lacerta  vivipara ),  respectively  (Lebreton  et  ah, 
1992),  and  time  has  been  used  to  model  the  survival  probability  of  the  northern  spotted 
owl  ( Strix  occidentalis  caurina)  (Forsman  et  ah,  1996)  and  the  North  Atlantic  right  whale 
(Caswell  et  ah,  1999). 

Covariates  are  incorporated  in  the  transition  probability  using  a  link  function.  The  link 
function  must  satisfy  the  constraint  that  each  column  of  the  transition  matrix  sums  to  1, 
and  each  entry  of  the  matrix  must  lie  between  0  and  1.  A  flexible  function  that  satisfies 
these  properties  is  the  polychotomous  logistic  function,  which  is  derived  by  expressing  the 
log  of  the  odds  ratio  as  a  linear  function  of  the  covariates  (Hosmer  &;  Lemeshow,  1989).  Let 
x[d^  be  the  value  of  dth  covariate  at  time  t.  The  polychotomous  logistic  function  is 


exp  (otjj  +  E dPj?^) 

1  +  Ej  exp  (otji  +  ZdPjfxtd)) 


(2.16) 


where  otji  is  an  intercept  parameter,  and  is  a  slope  parameter  associated  with  the 
dth  covariate.  When  all  the  slope  parameters  are  zero  for  all  d,  i,  and  j,  the  transition 
matrix  is  constant  over  time.  The  simple  logistic  function  that  is  often  used  in  mark- 
recapture  literatures  (e.g.  Burnham  et  ah,  1987;  Lebreton  et  ah,  1992)  is  a  special  case  of 
the  polychotomous  logistic  function  for  a  binary  outcome. 


2.4  Matrix  Population  Models 

Population  projection  matrices  contain  both  transition  probabilities  and  fertilities  (see 
Caswell,  2001).  Because  the  transition  probabilities  are  estimated  by  the  MSMR  method, 
we  can  construct  the  projection  matrix  if  we  know  the  fertility  terms.  In  this  section,  we 
show  an  example  of  how  those  terms  might  be  obtained,  and  how  to  compute  confidence 
intervals  for  population  growth  rate  calculated  from  the  population  matrix. 
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2.4.1  Conversion  from  a  Transition  Matrix  to  a  Population  Projection 


Matrix 


The  right  whale  example  provides  enough  information  to  write  a  2-sex  model.  To  do  so, 
we  renumber  the  male  stages  in  Figure  2.1  as  5,  6,  7.  Letting  4>ji(t)  denote  the  transition 
probability  as  before,  the  projection  matrix  is 


f 

0 

F2{t) 

F3(t) 

0 

0 

0 

0 

021  (0 

022(0 
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0 

0 

0 

032(0 

033(0 

034(0 

0 

0 

0 

0 

042  (0 

043(0 

0 

0 

0 

0 

0 

Fe(t) 

Fi{t) 

0 

0 

0 

0 

0 

0 

0 

0 

065(0 

066  (0 

0 

0 

0 

0 

0 

0 

076  it) 

077 (0  y 

(2.17) 


The  upper  left  and  lower  right  blocks  describe  production  of  females  by  females  and  males 
by  males,  respectively.  The  entries  in  the  lower  left  block  describe  production  of  males  by 
females. 

When  constructing  a  population  projection  matrix,  transition  probability  and  fertility 
terms  are  often  estimated  from  two  separate  data  sets  (Caswell,  2001),  but  the  fertility  terms 
can  be  directly  estimated  using  the  MSMR  method  if  the  stage  structure  includes  mothers 
that  give  birth  between  two  consecutive  sampling  periods  (i.e.  stage  4  in  our  example). 
Each  time  an  individual  enters  this  stage,  it  gives  birth;  therefore,  transition  probabilities 
into  the  fertile  stage  are  also  probabilities  of  giving  birth.  If  the  number  of  female  and 
male  births  at  each  reproductive  event  are  bf  and  bm,  respectively,  the  fertility  terms  in  the 
projection  matrix  are  given  by  the  product  of  the  number  of  offspring  and  the  transition 
probabilities: 


1*2  (0  =  6/042  (0> 


(2.18) 
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F3{t)  =  bf<t>4z(t), 

(2.19) 

Fe{t )  =  bm(j>42{t), 

(2.20) 

F7(t)  =  bm(f>43(t). 

(2.21) 

An  important  assumption  in  these  fertility  terms  is  that  mothers  and  their  newborns  have 
the  same  capturability  during  sampling.  To  ensure  that  this  condition  is  satisfied,  when 
mothers  are  captured,  their  offspring  should  also  be  captured  and  entered  into  the  data  as 
new  individuals.  Similarly,  when  newborns  were  captured,  their  mothers  should  be  captured 
and  identified  as  mothers.  Later,  we  will  show  one  example  of  remedial  methods  when  the 
equal  capturability  assumption  is  not  met. 

The  model  (2.17)  is  female-dominant;  males  do  not  affect  population  dynamics.  This 
assumption  is  often  legitimate  when  the  population  size  of  males  is  large  enough  that  search¬ 
ing  for  a  partner  does  not  limit  reproduction  by  females.  Thus  for  calculation  of  population 
growth  rate,  the  two-sex  matrix  may  be  reduced  to  the  female  matrix 


0 

m) 

Fs(t) 

0 

021  (t) 

022  (t) 

0 

0 

0 

032  (t) 

033(f) 

034(f) 

V  0 

042(f) 

043(f) 

0 

(2.22) 


2.4.2  Confidence  Intervals  for  Population  Growth  Rate 


The  long-term  population  growth  rate  implied  by  a  projection  matrix  At  is  given  by  the 
dominant  eigenvalue  A  of  Af.  A  confidence  interval  for  A  can  be  approximated  from  the 
MSMR  statistics,  using  eigenvalue  sensitivity  formula  and  the  covariance  matrix  of  the 
parameters. 

The  eigenvalue  sensitivity  formula  is 


<9  A  VjWi 
d  aji  <  w,  v  >  ’ 


(2.23) 
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where  v  and  w  are  the  left  and  right  dominant  eigenvectors  of  the  population  projection 
matrix  (Caswell,  1978,  2001).  If  the  a,ij  are  functions  of  some  other  parameters  0T,  the 
sensitivity  of  A  to  9r  is 

dX  i  dX  dcLjj  „  .n 

Now  let  9  be  a  vector  of  parameters  estimated  by  the  MSMR  method.  An  approximate 
95%  confidence  interval  for  A  is  calculated  by 


(2.25) 


where  cqr  is  (gr)th  entry  of  the  estimated  covariance  matrix  C.  The  covariance  matrix  C  can 
be  estimated  by  inverting  the  Hessian  matrix  (the  information  matrix)  (e.g.  Burnham  et  ah, 
1987).  This  method  of  constructing  the  confidence  interval  is  an  application  of  the  delta 
method  (see  Seber,  1982,  Chapter  1),  taking  advantage  of  the  existence  of  the  eigenvalue 
sensitivity  formula  (2.23). 


2.5  Application  to  the  North  Atlantic  Right  Whale 

We  have  applied  the  MSMR  method  to  data  on  the  North  Atlantic  right  whale  ( Eubalaena 
glacialis ).  The  northern  right  whale  is  considered  one  of  the  most  endangered  large  whale 
species  in  the  world  (Waring  et  al.,  1999).  The  current  population  in  the  western  North 
Atlantic  contains  fewer  than  300  individuals.  They  migrate  from  the  Bay  of  Fundy,  which  is 
a  summer  feeding  ground,  to  the  coast  of  Florida,  which  is  a  winter  calving  ground.  Caswell 
et  al.  (1999)  showed  that  the  crude  survival  probability  of  individuals  in  this  population 
has  been  declining  since  1980. 

Data  on  the  North  Atlantic  right  whale  have  been  collected  by  New  England  Aquarium 
and  consist  of  annual  sighting  histories  of  photographically  identified  animals  from  1980  to 
1997  (Crone  &  Kraus,  1990).  For  the  purpose  of  our  analysis,  we  consider  individuals  to 
have  been  marked  on  the  occasion  of  their  first  identification,  and  recaptured  when  they 
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were  resighted  during  a  subsequent  year.  Of  the  372  individuals  used  for  the  analysis,  141 
are  known  to  be  females  and  143  to  be  males.  We  assumed  the  remainder  to  be  either 
female  or  male  with  50%  probabilities.  A  few  sightings  of  dead  individuals  exist  but  are 
not  included  in  this  analysis. 

2.5.1  Stage-assignment  matrices 

We  attempted  to  assign  each  individual  at  each  capture  to  one  of  the  stages  shown  in 
Fig.  2.1.  A  whale  was  considered  mature  if  it  was  known  to  be  at  least  9  years  old  or, 
for  females,  if  it  had  been  observed  with  a  calf.  Stages  could  be  assigned  with  certainty 
in  78%  of  the  captures.  The  remainder  were  known  to  be  either  immature  or  mature,  and 
for  these  captures,  we  must  calculate  the  entries  U2,  U3,  u$,  and  uq  of  the  stage-assignment 
matrices  (2.5)  and  (2.6).  In  the  absence  of  information  to  the  contrary,  we  assume  that  these 
probabilities  are  constant  over  time  and  across  individuals,  but  differ  between  females  and 
males.  Because  we  use  different  criteria  to  assign  females  and  males  to  stages,  we  expect 
that  the  probability  distribution  of  stages  among  the  unknown  staged  captures  would  differ 
for  females  and  males. 

When  the  stage  of  a  captured  individual  is  uncertain,  the  (2,2)  entry  of  the  stage- 
assignment  matrix  is  the  probability  that  the  individual  is  immature,  given  that  the  stage 
is  uncertain.  Similarly,  the  (3,3)  entry  is  the  probability  that  the  individual  is  mature,  given 
that  the  stage  is  uncertain.  The  other  entries  are  all  zero.  To  express  these  probabilities  in 
mathematical  form,  let  X  be  a  random  variable  giving  the  stage  of  an  individual  and  Y  be 
a  random  variable  taking  the  value  1  if  the  stage  is  known  and  0  if  the  stage  is  uncertain. 
Then  the  two  probabilities  are 


Pr{X  =  2|  Y  =  0) 

(2.26) 

Pr{X  =  3\Y  =  0) 

1  -  U2- 

(2.27) 
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To  calculate  Pr(X  =  2| Y  =  0),  we  use  Bayes’  Rule  to  derive: 


Pr(X  =  2\Y  =  0)  = 


Pr(X  =  2)  -  Pr{X  =  2|  Y  =  1  )Pr(Y  =  1) 
1  -  Pr(Y  =  1) 


(2.28) 


Here,  Pr(Y  =  1)  is  the  probability  that  the  stage  of  a  captured  immature  or  mature 
individual  is  known,  and  can  be  estimated  from  the  capture  history  data  as 


Pr(Y  =  1)  = 


N2  +  N3 
Ni  +  Nz  +  Nu’ 


(2.29) 


where  N2,  N3,  and  Nu  are  numbers  of  captures  of  immature,  mature,  and  uncertain  stages, 
respectively.  Pr(X  =  2| Y  =  1)  is  the  probability  that  the  stage  of  a  detected  immature 
or  mature  individual  is  immature  given  that  the  stage  is  known.  This  probability  can  be 
calculated  from  the  capture  history  data  as 

Pr(X  =  2|  Y  =  1)  =  (2.30) 


Finally,  Pr(X  =  2)  is  the  probability  that  the  stage  is  immature  given  that  the  stage  is  either 
immature  or  mature,  regardless  of  whether  the  stage  is  known  or  uncertain.  To  estimate 
this  probability,  we  estimated  the  parameters  for  a  time-invariant  projection  matrix  from 
the  subset  of  the  data  containing  only  certain  captures.  From  the  stable  stage  distribution 
w  (i.e.  the  right  eigenvector  associated  with  the  dominant  eigenvalue)  of  this  matrix,  we 
calculated  the  proportion  of  individual  in  stage  2  among  stages  2  and  3  and  used  it  as  our 
estimate  of  Pr(X  =  2): 


Pr(X  =  2)  = 


w  2 

W2  +  U>3 


(2.31) 


For  males,  the  same  method  was  applied  to  the  male  stages  (5  and  6).  It  should  be  noted 
that  these  calculations  work  best  when  the  capture  probability  of  stages  2  and  4  (5  and  6 
for  males)  are  similar.  Otherwise,  each  count  in  (2.29)  and  (2.30)  should  be  divided  by  the 
corresponding  capture  probability  (Nichols  et  al.,  1994).  The  end  result  of  these  calculations 
is  in  U2  =  0.87,  uz  =  0.13,  1x5  =  0.30,  and  ue  —  0.70. 
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2.5.2  Capture  probabilities 

Capture  probabilities  were  modeled  as  binary  logistic  functions  of  estimated  sampling  effort 
levels  in  the  Northern  and  Southern  regions,  which  are  major  feeding  and  calving  grounds, 
respectively.  These  effort  levels  were  approximated  by  the  number  of  sampling  dates  per  year 
in  each  region.  We  created  models  by  including  all  possible  combinations  of  effort  levels  for 
all  possible  combination  of  stages.  This  resulted  in  1024  models.  The  best  capture  model 
among  the  1024  candidate  models  was  selected  using  Akaike  Information  Criteria  (AIC; 
Akaike  (1973)).  Because  the  sample  size  is  large,  we  did  not  use  the  small  sample  adjustment 
to  AIC  (i.e.  AICc  in  Burnham  &  Anderson  (1998)).  The  differences  in  AIC  between  the 
best  and  the  second-best  capture  models  was  about  2,  indicating  that  the  support  for 
the  best  model  relative  to  the  second  best  model  is  high  (Burnham  &  Anderson,  1998). 
Furthermore,  the  4  best  models  differ  only  in  how  capture  probability  of  mothers  depends 
on  effort,  and  in  all  cases  the  capture  probability  was  consistently  close  to  1  throughout  the 
sampling  period.  Therefore,  we  used  only  the  best  model  shown  in  Table  2.  The  capture 
probabilities  of  immature  males  and  females  did  not  differ  significantly  in  the  best  model 
based  on  a  likelihood  ratio  test.  Therefore,  we  set  these  two  capture  probabilities  equal  and 
used  the  resulting  capture  probability  model  for  further  analysis  (Table  2.5.2). 


Table  2.2:  Dependence  of  the  best  capture  model  for  the  North  Atlantic  right  whale  on 
effort  level  and  time. 


Stage* 

Northern  Effort 

Southern  Effort 

Immature  Female  (2) 

Yes 

Yes 

Mature  Female  (3) 

Yes 

Yes 

Mature  Female  with  Calf  (4) 

No 

No 

Immature  Male  (6) 

Yes 

Yes 

Mature  Male  (7) 

Yes 

Yes 

*  Sighting  probability  of  calves  cannot  be  estimated  because  the  capture 
of  a  calf  is  always  the  first  capture  of  that  individual. 
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2.5.3  Transition  Probabilities 


Although  we  have  shown  the  vital  rates  have  varied  over  time  (Chapter  5),  for  this  example 
we  fit  a  model  in  which  the  transition  probabilities  are  constant  over  time  (i.e.  no  covariates). 
We  also  assumed  that  the  survival  probabilities  of  female  and  male  calves  are  the  same. 
This  model  gives  a  time-averaged  picture  of  right  whale  demography.  Estimated  capture 
and  transition  probabilities  are  shown  in  Fig.  2.2  and  Table  2.3,  respectively. 

The  population  projection  matrix  for  female  right  whales  is 


f  0  0.5^42-^34  0.5^431/^34  0 

<t>2\  0  22  0  0 

0  0  32  </>33  0  34 

^  0  0 42  043  0  ) 


(2.32) 


Table  2.3:  Estimated  transition  probabilities  for  the  North  Atlantic  right  whale. 


Transition 

Probability 

Mean 

95%  Confidence 

Interval 

021 

0.92 

[0.74  0.98] 

022 

0.86 

[0.81  0.89] 

032 

0.08 

[0.06  0.12] 

042 

0.02 

[0.01  0.03] 

033 

0.80 

[0.77  0.83] 

043 

0.19 

[0.16  0.22] 

034 

0.83 

[0.77  0.88] 

066 

0.76 

[0.72  0.79] 

076 

0.19 

[0.16  0.23] 

077 

0.95 

[0.94  0.96] 

Note:  The  confidence  intervals  were  estimated 

from  1000  parametric  bootstrap  samples  generated 

assuming  multivariate  normal  distributions  of 

parameters.  For  the  covariance  matrix  of  the 

distribution  was  estimated  as  the  inverse  of 

the  Hessian  matrix  (Burnham  et  al.,  1987;  Lebreton,  1995). 
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Year 


Figure  2.2:  Stage-specific  capture  probabilities  for  (a)  immature  male  and  female,  (b)  ma¬ 
ture  female,  and  (c)  mature  male  right  whales.  Error  bars  indicate  point-wise  95%  con¬ 
fidence  intervals  estimated  from  1000  parametric  bootstrap  samples  generated  assuming  a 
multivariate  normal  distribution  of  the  logit  of  parameters.  The  covariance  matrix  of  the 
distribution  was  estimated  as  the  inverse  of  the  Hessian  matrix  (see  Burnham  et  al.,  1987; 
Lebreton,  1995).  Mothers  had  a  constant  capture  probability  0.99  (95%  Cl  =[0.98,  1.00]). 


35 


This  matrix  is  the  same  as  (2.22)  but  with  a  particular  set  of  assumptions  defining  the 
fertility  terms.  Consider  i^t)  in  (2.22).  When  a  female  moves  from  stage  2  to  stage  4 
(with  probability  <^42),  she  gives  birth;  the  newborn  is  female  with  probability  0.5.  To 
appear  as  a  calf  in  stage  1  at  t  +  1,  the  newborn  calf  must  survive  long  enough  to  be 
cataloged.  Although  newborn  calves  have  distinct  markings,  they  are  harder  to  distinguish 
individually  than  other  stages.  Therefore,  calf  survival  is  estimated  from  the  time  when 
the  calf  is  seen  sufficiently  well  to  permit  identification,  which  is  not  necessarily  on  its  first 
sighting.  We  assumed  calves  are  identified  on  average  midway  through  their  first  year, 
and  that  the  mother  must  survive  this  long  (with  probability  )  in  order  for  the  calf  to 
survive.  Fs(t)  is  calculated  in  a  similar  manner. 

Prom  this  matrix,  we  estimated  the  long  term  population  growth  rate  and  its  confidence 
interval.  They  are  A  =  1.01  (95%  Cl  =[1.00,  1.02]).  This  result  shows  that  the  North 
Atlantic  right  whale  population  has  been  growing  by  1%  annually  on  average  from  1980  to 
1997.  (In  fact,  a  time- varying  model  estimated  by  this  same  procedure  concludes  that  the 
growth  rate  has  declined  from  A  =  1.03  to  A  =  0.98  over  this  time  period  Chapter  5)  This 
matrix  can  now  be  analyzed  to  obtain  the  stable  stage  distribution,  reproductive  value, 
damping  ratio,  sensitivity  and  elasticity  of  A,  and  other  demographic  statistics. 

2.6  Discussion 

The  method  presented  here  estimates  a  population  projection  matrix  from  mark-recapture 
data,  which  is  one  of  the  most  commonly  available  data  types  for  animal  populations. 
Once  the  population  projection  matrix  is  estimated,  it  is  subject  to  complete  demographic 
analysis;  such  analyses  provide  powerful  tools  for  conservation  biology  (e.g  Caswell,  1989; 
Tuljapurkar  &  Caswell,  1997;  Caswell,  2001).  They  can  be  used  to  assess  the  causes  of 
past  population  declines,  and  to  predict  the  effect  of  possible  future  management  actions. 
Because  population  projection  matrices  contain  many  parameters,  it  has  been  difficult  to 
estimate  them  accurately.  This  has  been  especially  true  for  animals  that  are  not  captured 
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at  every  sampling  period. 

The  likelihood  calculations  here  are  simpler  than  those  described  in  Nichols  et  al. 
(1992b).  This  allows  the  use  of  mathematical  software  packages  such  as  Matlab,  so  the 
transition  and  capture  probability  models  need  not  be  limited  to  those  available  in  mark- 
recapture  packages  such  as  Mark  (White  &  Burnham,  1999)  ,  Mssurviv  (Hines,  1994),  or 
Surviv  (White,  1983). 

Our  method  permits  the  use  of  capture  histories  with  uncertain  stage  and  sex.  Individ¬ 
uals  with  such  uncertainties  tend  to  have  lower  survival  rates  than  the  rest  of  a  population 
because  individuals  that  survive  longer  have  more  chances  for  accurate  assessment  of  stage 
and  sex  identification.  For  example,  right  whales  are  considered  mature  at  9  years  of  age. 
If  animals  that  die  within  9  years  from  their  first  capture  are  excluded  because  their  stage 
is  uncertain,  then  we  would  overestimate  the  survival  probability.  Observations  with  un¬ 
certain  stage  or  sex  should  not  be  discarded  in  parameter  estimation.  Our  approach  is  one 
way  to  deal  with  this  problem. 

The  stage  structure  we  used  in  this  paper  contains  as  a  stage  females  that  have  given 
birth  between  consecutive  sampling  periods.  This  stage  makes  the  conversion  of  the  tran¬ 
sition  matrix  into  a  population  projection  matrix  relatively  simple.  Because  the  purpose  of 
the  MSMR  statistics  often  is  to  estimate  a  population  projection  matrix,  we  recommend 
the  use  of  this  type  of  stage  structure  when  possible. 

The  polychotomous  logistic  function  is  a  flexible  way  to  allow  transition  probabilities 
to  decrease  or  increase  with  a  covariate  while  satisfying  the  requirement  that  each  column 
of  the  transition  matrix  sum  to  one.  When  time  is  used  as  a  covariate,  the  polychotomous 
function  allows  inferences  about  temporal  trends  in  stage-specific  transition  rates.  This 
approach  has  been  applied  to  the  North  Atlantic  right  whale  data  (Chapter  5). 

Multi-stage  mark-recapture  data  arise  in  many  applications.  For  example,  Nichols  & 
Kendall  (1995)  use  them  in  population  genetics  context  to  test  trade-offs  between  survival 
and  reproduction.  Hestbeck  et  al.  (1991)  use  them  to  estimate  spatial  movement  of  indi¬ 
viduals.  We  have  applied  them  to  deal  with  the  problem  of  temporary  emigration  (Chapter 
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3).  We  hope  the  extentions  of  the  analytical  method  presented  here  will  make  them  even 
more  useful. 

Mark-recapture  data  are  expensive  to  collect,  and  they  should  be  analyzed  as  completely 
as  possible.  If  information  on  the  stage  of  individuals  (e.g.,  age,  size,  other  developmental 
stages,  or  geographic  locations)  is  collected  in  addition  to  the  basic  mark-recapture  data, 
then  MSMR  statistics  can  be  applied.  The  stage  information  need  not  be  complete  because 
our  method  incorporates  uncertainties  in  stage  identifications.  The  value  of  being  able  to 
use  matrix  population  models  for  conservation  makes  it  worthwhile  to  collect  stage-specific 
data  whenever  possible. 
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Chapter  3 


A  general  approach  to  temporary 
emigration  in  mark-recapture 
analysis 

3.1  Introduction 

In  mark-recapture  studies,  survival  probability  is  estimated  from  capture  histories  of  in¬ 
dividually  identified  animals.  Those  histories  contain  information  on  whether  or  not  each 
individual  was  captured  at  each  sampling  occasion.  The  analysis  assumes  that  all  individu¬ 
als  have  identical  capture  probabilities.  This  assumption  is  violated  when  some  individuals 
temporarily  leave  the  sampling  area  and  return  during  subsequent  sampling  occasions.  We 
refer  to  this  as  “temporary  emigration.”  The  capture  probability  of  the  individuals  that 
have  emigrated  is  zero,  while  the  rest  of  individuals  that  are  alive  have  a  non-zero  capture 
probability.  Therefore,  assuming  that  all  individuals  have  the  same  capture  probability, 
regardless  their  location,  results  in  under-  and  over-estimation  of  the  capture  probabil¬ 
ity  for  individuals  inside  and  outside  the  sampling  area,  respectively.  The  biased  capture 
probability  estimate  often  biases  survival  probability  estimates. 

The  temporary  emigration  process  we  consider  in  this  paper  is  deterministic;  i.e.,  we 
assume  that  all  animals  in  the  sampling  area  emigrate,  with  probability  1,  before  the  next 
sampling  occasion.  This  type  of  temporary  emigration  process  is  common  during  imma¬ 
ture  stages  and  between  reproductive  events.  We  refer  to  these  as  the  “immature  emigra- 
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tion  process”  and  “inter-birth  emigration  process,”  respectively.  For  example,  albatrosses 
(Diomedeidae  spp.)  leave  their  breeding  islands  immediately  after  fledging  (immature  em¬ 
igration  process)  and  do  not  return  for  5-15  years.  Once  they  return,  they  are  recaptured 
only  when  they  breed  (inter-birth  emigration  process)  (Weimerskirch  et  al.,  1997).  Simi¬ 
larly,  grey  seals  ( Halichoerus  gryprus )  are  marked  when  young,  but  cannot  be  recaptured 
before  they  mature  (immature  emigration  process)  (Schwarz  &  Stobo,  2000).  Right  whales 
( Eubalaena  sp)  are  monitored  at  their  calving  ground,  which  they  leave  after  reproduction 
and  return  at  intervals  of  several  years  for  calving  (inter-birth  emigration  process)  (Payne 
et  al.,  1990).  Similarly,  sea  turtles  are  often  individually  identified  only  at  their  natal  beach 
when  they  come  back  to  lay  their  eggs  (inter-birth  emigration  process). 

Deterministic  temporary  emigration  is  different  from  permanent  emigration  and  purely 
random  temporary  emigration,  which  are  also  common.  Permanent  emigration  occurs  when 
animals  leave  the  sampling  area  and  never  come  back.  In  usual  mark-recapture  analyses, 
permanent  emigration  is  either  assumed  to  be  a  part  of  mortality  as  a  loss  from  a  population 
(e.g.  Lebreton  et  al.,  1992)  or  assumed  not  to  occur  (e.g.  Caswell  et  al.,  1999).  In  purely 
random  temporary  emigration,  individuals  leave  and  return  independently  of  their  location 
at  the  previous  sampling  occasion.  This  type  of  temporary  emigration  occurs  commonly 
when  only  a  part  of  a  habitat  is  sampled.  When  mark- recapture  methods  are  applied 
to  a  population  exhibiting  purely  random  temporary  emigration,  the  estimated  capture 
probability  is  the  product  of  the  probability  that  individuals  are  in  the  sampling  area  and 
the  probability  that  individuals  within  the  sampling  area  are  captured.  However,  survival 
probability  is  not  affected  by  the  emigration  process  (Brownie  et  al.,  1993). 

In  this  paper,  we  show  how  to  estimate  survival  probability  and  other  demographic 
transitions  in  the  presence  of  deterministic  temporary  emigration.  Our  method  includes 
location  (temporarily  emigrated  or  not)  as  a  property  of  the  individual  and  uses  a  multi¬ 
stage  mark-recapture  method.  Our  approach  is  related  to  the  methods  used  by  Arnason 
(1972),  Arnason  (1973),  and  Hestbeck  et  al.  (1991)  to  estimate  movement  rates,  but  applies 
to  populations  that  exhibit  temporary  emigration.  Because  individuals  that  have  emigrated 


are  unobservable,  the  capture  history  data  usually  do  not  contain  enough  information  to 
estimate  all  parameters  separately.  In  order  to  add  necessary  information,  we  introduce  con¬ 
straints  among  the  parameters.  These  constraints  may  be  based  on  biologically  reasonable 
models  (e.g.,  that  temporary  emigration  does  not  influence  survival),  or  may  be  provided  by 
independent  estimates  of  some  parameters  from  other  data  sets.  The  constraints  reduce  the 
number  of  parameters  to  be  estimated,  and  allow  us  to  estimate  all  the  remaining  parameter 
values  separately.  The  possibility  of  modeling  mark-recapture  data  with  stage  structures 
that  contain  unobservable  stages  was  suggested  by  Lebreton  (1995).  Here,  we  apply  the 
idea  to  populations  that  exhibit  the  two  types  of  temporary  emigration  and  analyze  the 
estimability  of  parameters.  Our  results  generalize  those  of  Clobert  et  al.  (1994),  Lebreton 
et  al.  (1999),  and  Pradel  &  Lebreton  (1999)  in  that  we  apply  the  method  to  inter-birth 
temporary  emigration,  analyze  more  general  cases  of  immature  temporary  emigration,  and 
use  a  formal  method  to  determine  estimable  parameters. 

As  usual,  estimating  stage-specific  transition  and  capture  probabilities  involves  con¬ 
structing  an  appropriate  stage  structure,  expressing  the  likelihood  function  in  terms  of 
the  parameters,  and  finding  the  best  fit  parameter  values  using  maximum  likelihood  (e.g. 
Nichols  et  ah,  1992b,  also  Chapter  2).  The  likelihood  function  can  be  written  down  directly 
using  Markov  chains  that  specify  transition  probabilities  among  stages,  assuming  indepen¬ 
dence  of  fates  among  individuals  (i.e.  assuming  the  number  of  outcomes  falling  into  the 
possible  capture  history  sequences  is  multinomial)  (Caswell,  2001,  also  Chapter  2).  This 
likelihood  can  be  maximized  numerically  to  find  the  best  fit  parameter  values.  This  method 
allows  simple  programming  using  software  such  as  Matlab.  Alternatively,  in  many  but 
not  all  cases,  parameters  can  be  estimated  using  softwares  that  are  specifically  designed  for 
mark-recapture  analyses  such  as  Mark  (White  &;  Burnham,  1999),  MSSURVIV  (Hines, 
1994),  and  Surviv  (White,  1983). 
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3.2  Temporary-emigration  stage  structures 

Figures  3.1  and  3.2  show  examples  of  life  cycle  graphs  including  temporary  emigration. 

Each  arrow  indicates  a  possible  stage  transition  from  time  t  to  t  + 1  and  has  a  transition 
probability  associated  with  it.  We  denote  the  transition  probability  from  stage  i  to  j  by 
<j)ji.  It  should  be  noted  that  the  subscripts  for  the  transition  probability  are  reversed  from 
the  traditional  notation  (e.g.  Amason,  1973;  Hestbeck  et  al.,  1991)  so  that  they  indicate  the 
entry  in  a  matrix  population  model  (e.g.  Caswell,  2001)  and  the  column  stochastic  transition 
matrix  used  in  Chapter  2.  Individuals  in  a  capturable  stage  have  a  non- zero  probability  of 
being  recaptured.  This  probability  is  denoted  by  pi  for  stage  i. 

The  stage  structure  for  the  inter-birth  emigration  process  (Fig.  3.1a)  includes  a  breeding 
adult  stage  (stage  1)  in  which  individuals  can  be  recaptured  with  a  probability  pi,  and 
k  —  1  inter-breeding  adult  stages  in  which  individuals  cannot  be  recaptured.  The  value  of 
k  depends  on  the  minimum  inter-birth  interval.  For  example,  Wandering  Albatrosses  never 
breed  successfully  in  two  consecutive  years  (Weimerskirch  et  al.,  1997).  Those  that  give 
birth  in  one  year  will  be  missing  from  the  breeding  ground  in  the  following  year.  After  two 
years,  they  return  to  the  breeding  ground  with  some  probability  that  is  independent  of  age. 
Therefore,  k  =  2  for  Wandering  Albatross  (Fig.  3.1b).  Similarly,  individual  right  whales  do 
not  give  birth  for  at  least  two  years  after  successful  reproduction,  so  k  =  3  for  the  right 
whale  (Fig.  3.1c). 

The  stage  structure  for  the  immature  emigration  process  (Fig.  3.2a)  contains  a  newborn 
stage  (stage  0)  in  which  individuals  are  first  captured  and  marked,  k  —  1  immature  stages 
(stage  1  through  k- 1)  in  which  individuals  cannot  be  captured  due  to  temporary  emigration, 
and  a  mature  stage  (stage  k )  in  which  individuals  can  be  captured  with  a  probability  pk- 
If  k  >  1,  the  first  k  -  1  stages  (i.e.  newborn  and  the  first  k-  2  immature  stages)  represent 
individuals  that  cannot  mature  before  the  following  sampling  period,  and  these  stages  also 
correspond  to  ages  of  individuals.  Stage  A:  —  1  may  mature  (with  probability  <pk,k- 1),  remain 
immature  (with  probability  i,fc— l),  or  die  (with  probability  1  —  {<j>k-\,k-i  +  <Pk,k- 1))- 
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Figure  3.1:  Inter-birth  temporary-emigration  stage  structures,  (a)  General  stage  structure, 
(b)  Stage  structure  with  k  =  1.  This  represents  the  temporary  emigration  of  Wandering 
Albatross.  (Models  1.1,  1.2,  and  1.3).  (b)  Stage  structure  with  k  =  2.  This  represents 
the  temporary  emigration  of  the  right  whale.  (Models  2.1,  2.2,  and  2.3).  Parameters  for 
these  stage  structures  are  fcj  a  transition  probability  from  stage  j  to  i,  and  pi  a  capture 
probability  of  stage  1.  Stages  other  than  stage  1  have  a  zero  recpature  probability. 
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Figure  3.2:  Immature  temporary-emigration  stage  structure,  (a)  General  stage  structure, 
(b)  Stage  structure  with  k  =  4.  This  represents  the  temporary  emigration  of  the  grey  seal. 
(Models  3.1,  3.2,  3.3,  3.4,  and  3.5).  (c)  Stage  structure  with  k  =  5.  This  represents  the 
temporary  emigration  of  Wandering  Albatross.  (Models  4.1,  4.2,  4.3,  and  4.4).  Parameters 
for  these  stage  structures  are  <j>ij  the  transition  probability  from  stage  j  to  i,  and  pk  the 
recapture  probability  of  the  mature  stage.  Stages  other  than  stage  k  have  a  zero  recpature 
probability. 
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For  example,  grey  seals  do  not  give  birth  before  reaching  age  4  (Schwarz  &  Stobo,  2000). 
Therefore,  an  appropriate  stage  structure  for  them  is  one  with  k  =  4  (Fig.  3.2b).  Similarly, 
albatrosses  do  not  mature  before  reaching  the  age  of  5  (Weimerskirch  et  ah,  1997);  therefore, 
k  —  5  in  their  stage  structure  (Fig.  3.2c).  When  k  =  1,  the  stage  structure  does  not  contain 
a  separate  immature  stage  (i.e.  newborn  and  immature  stages  are  the  same),  and  newborn 
individuals  can  mature  directly  or  remain  immature.  In  this  stage-structure,  immature 
individuals  can  be  marked  but  cannot  be  recaptured;  therefore,  they  still  exhibit  temporary 
emigration.  These  stage  structures  include  as  special  cases  the  model  of  Clobert  et  al. 
(1994),  but  relax  their  assumption  of  a  maximum  age  in  immature  stage. 

Other  population  parameters  (e.g.,  stage- specific  survival  probabilities,  stage-specific 
probabilities  of  death,  and  probabilities  of  breeding)  can  be  calculated  from  the  faj .  For 
examples,  the  stage  specific  survival  probability  (s{)  is  given  by  4>ji  and  the  stage  specific 
probability  of  death  is  given  by  1  -  <t> ji  ■  In  the  immature  emigration  model,  the  breeding 
probability,  conditional  on  survival  of  individuals  in  stage  k  —  1  (we  denote  this  probability 


by  VO  is 


V>  = 


<f>k,k- 1 


<Pk,k- 1  +  4>k- 1, 

We  will  show  how  to  calculate  some  other  population 
be  found  in  (Caswell,  2001). 


— .  (3.1) 

k- 1 

indices  in  Section  6.  More  details  can 


3.3  Estimability  of  parameters 

Unfortunately,  some  of  the  parameters  in  the  stage  structures  for  temporary  emigration 
processes  may  not  be  estimable,  because  of  the  lack  of  observations  of  the  emigrated  stages. 
In  general,  whether  best-fit  parameter  values  can  be  found  uniquely  or  not  depends  on  how 
the  parameters  appear  in  the  likelihood  function.  If  two  or  more  parameters  are  confounded 
(i.e.  they  appear  only  in  the  same  arithmetic  form  throughout  the  likelihood  function),  they 
cannot  be  estimated  separately.  The  solution  is  to  impose  constraints  on  the  parameters 
by  specifying  a  model  for  the  relationships  among  them. 
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Confounded  parameters  may  result  from  model  specification  (intrinsic  identifiability 
problem)  or  a  lack  of  variation  among  observed  data  (extrinsic  identifiability  problem) 
(McCullagh  &  Nelder,  1989).  When  parameters  are  intrinsically  confounded,  the  only  way 
to  eliminate  the  problem  is  to  constrain  parameters  or  modify  the  stage-structure  so  that 
other  useful  parameters  can  be  estimated  separately.  When  parameters  are  confounded 
because  of  lack  of  variation  among  observed  data,  the  solution  is  to  constrain  parameters, 
modify  the  stage-structure,  or  increase  the  sample  size.  We  will  show  how  to  determine 
whether  model  has  the  intrinsic  and  extrinsic  identifiability  problems  or  not  in  the  next 
section  and  give  examples  in  Section  4  and  5. 

3.3.1  Determining  estimability  of  parameters 

Whether  parameters  are  confounded  in  the  likelihood  function  or  not  can  be  determined 
using  the  method  of  Catchpole  &  Morgan  (1997).  Their  method  uses  the  Jacobian  matrix, 
which  contains  the  derivatives  of  the  likelihood  functions  with  respect  to  parameters.  If  the 
rank  of  the  Jacobian  equals  the  number  of  parameters,  then  all  parameters  can  be  estimated 
separately.  Here,  we  give  a  heuristic  description  of  the  method  as  it  applies  to  temporary 
emigration  models;  for  details  see  Catchpole  &  Morgan  (1997). 

We  start  with  data  consisting  of  capture  histories  of  n  individuals.  The  contribution 
of  individual  i  to  the  likelihood  function  is  proportional  to  the  probability  of  its  capture 
history.  Let  k  be  this  contribution  and  assume  that  the  model  contains  m  parameters  (aq, 
X2,  ...,  Xm).  Then,  k  is  a  function  of  the  parameters: 

li  =  X2,  Xnx).  (3.2) 

Taking  the  partial  derivatives  of  the  function  k  with  respect  to  the  parameters  results  in 

d  li  =  Tr^-dxi  +  dx2  +  ...  +  —  dxm.  (3-3) 

OX\  UX  2  £'*•£771 

Because  the  data  consist  of  n  captured  individuals,  n  such  functions  exist.  In  a  vector 
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notation,  they  are: 


dl  =  Jdx 


(3.4) 


where, 


d,=l 

dl\  dh  dfa 

•  •  dln 

dx  = 

^  dxi 

dX2  •  • 

dxm 

(  dR 

dx\ 

dh. 

dx2 

dh  \ 

dxm 

dj± 

dx\ 

dh 

8x2 

dh 

dx-m 

J  = 

dh 

dxi 

dh 

dx  2 

dh 

dx-m 

\  it 

dfn 

dl2 

OXm  / 

Here,  J  is  the  Jacobian  matrix  of  the  likelihood  function  with  respect  to  the  parameters. 
Equation  (3.4)  is  a  system  of  linear  equations  with  J  as  a  coefficient  matrix.  For  x  to 
be  estimated  uniquely,  the  Jacobian  matrix  must  map  dx  to  dl  uniquely.  This  requires 
that  the  rank  of  the  Jacobian  matrix  be  m  (i.e.  columns  of  the  Jacobian  matrix  must  be 
independent).  Therefore,  if  the  rank  of  J  is  equal  to  the  number  of  parameters,  values  for  the 
parameters  can  be  estimated  separately  using  the  maximum  likelihood  method.  Otherwise, 
some  parameters  cannot  be  estimated  separately,  and  the  number  of  parameters  that  can 
be  estimated  separately  is  given  by  the  rank  of  J. 

If  multiple  individuals  have  the  same  capture  history  sequence,  the  likelihoods  associated 
with  their  capture  histories  are  the  same.  Therefore,  only  one  such  copy  needs  to  be 
included  in  the  likelihood  vector,  and  the  duplicated  likelihoods  can  be  eliminated.  If  we 
do  the  analysis  assuming  that  every  possible  capture  history  is  observed  at  least  once  for 
each  stage  structure,  then  the  above  analysis  can  be  used  to  detect  intrinsic  identifiability 
problems.  On  the  other  hand,  if  we  use  likelihood  associated  with  the  available  capture 
history  sequences,  the  analysis  can  be  used  to  detect  extrinsic  identifiability  problems. 
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3.3.2  An  example  of  the  Jacobian  method 

As  an  example,  we  apply  the  above  method  to  the  inter-birth  emigration  model  with  k  =  2 
(Fig.  3.1b).  We  assume  all  individuals  are  marked  at  the  first  sampling  occasion,  and  there 
are  4  subsequent  sampling  occasions  for  recaptures  (i.e.  5  capture  occasions).  Under  these 
assumptions,  there  are  five  possible  capture  histories: 

hi  =  10101, 

h2  =  10100, 
hz  =  10001, 
h4  =  10010, 
h5  =  10000, 

where  1  indicates  an  individual  was  captured  and  0  indicates  it  was  not  captured. 

The  likelihood  functions  corresponding  to  these  capture  histories  are: 


h  = 

Pi  021 012  > 

(3.5) 

h  — 

Pi  021012  —  Pi  021012 

(3.6) 

h  = 

[021012(1  -  pi)  +  021022]012Pl 

(3.7) 

h  = 

Pi  021 022012  [021  +  (1  —  021 )] 

(3.8) 

h  = 

1  —  h  —  l2  —  h  —  h- 

(3.9) 

A  likelihood  function  can  be  determined  by  creating  a  list  of  event  sequences  that  could 
lead  to  the  capture  history,  calculating  the  probability  of  the  sequences,  and  then  summing 
the  probabilities.  For  example,  hi  can  only  arise  for  an  animal  that  makes  the  transition 
back  to  stage  1  immediately  after  each  observation  occasion  in  stage  2,  and  is  observed 
on  each  occasion  it  is  in  stage  1.  The  corresponding  likelihood  function  (equation  3.5) 
is  simply  the  product  of  this  unique  sequence  of  transition  and  recapture  probabilities, 
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with  some  collection  of  terms.  Likelihood  functions  1 2  through  I4  are  more  complicated 
because  more  than  one  sequence  of  events  can  lead  to  the  observation  sequence  /12  through 
hi.  For  example,  /13  could  arise  from  an  animal  that  remained  in  stage  2  during  sampling 
occasions  2-4  and  then  made  the  transition  to  stage  2  for  sampling  occasion  5  and  was 
recaptured.  However,  this  sequence  could  also  arise  if  the  animal  went  through  the  same 
stage  transitions  that  produced  hi,  but  was  not  recaptured  during  sampling  occasion  3.  l5 
is  calculated  simply  as  the  probability  that  capture  histories  hi  through  /14  do  not  occur. 

There  are  4  parameters  (pi,  (f> 21,  (p22,  and  <£>12)  and  5  likelihoods;  therefore,  the  Jacobian 
is  a  5x4  matrix.  Its  rank  is  the  number  of  independent  columns  in  J.  The  Jacobian  and  its 
rank  can  be  found  easily  using  software  for  symbolic  calculations;  we  used  the  Symbolic 
Math  Toolbox  in  Matlab.  The  rank  of  the  Jacobian  is  3,  so  only  3  of  the  4  parameters 
can  be  estimated  separately;  therefore,  this  model  has  an  intrinsic  identifiability  problem. 
This  does  not  necessarily  imply  that  any  3  parameters  can  be  estimated;  we  include  a 
discussion  on  which  3  parameters  can  be  estimated  in  the  next  section. 

3.4  Estimable  parameters  in  temporary  emigration  models 

To  evaluate  estimability  of  parameters  in  selected  temporary  emigration  models,  we  applied 
the  Jacobian  method  to  15  models  based  on  the  four  stage  structures  in  Fig.  3.1b,  3.1c,  3.2b, 
and  3.2c.  Table  3.1  lists  the  models  and  also  shows  (1)  constraints  on  transition  probabilities, 
(2)  parameters  whose  values  are  assumed  to  be  known,  (3)  the  number  of  parameter  to 
be  estimated,  and  (4)  the  rank  of  the  Jacobian  matrix.  After  imposing  constraints  and 
eliminating  known  parameters,  all  the  remaining  parameters  become  estimable  in  some 
models,  which  are  indicated  by  (*). 

Consider  the  inter-birth  emigration  model  of  Wandering  Albatross  (Fig.  3.1b,  k  =  2). 
Adults  can  be  captured  on  the  breeding  grounds,  but  then  disappear  until  they  breed  next. 
If  no  constraints  are  placed  on  the  parameters  (Model  1.1),  the  rank  of  the  Jacobian  is  3, 
so  we  can  estimate  only  at  most  3  of  the  4  parameters.  However,  this  does  not  imply  that 
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any  3  parameters  can  be  estimated.  For  example,  if  pi  is  known  independently  (Model  1.2), 
the  number  of  parameters  is  reduced  to  3,  but  the  rank  of  the  Jacobian  is  reduced  to  2. 
On  the  other  hand,  we  can  fit  a  model  that  assumes  that  survival  probability  of  breeding 
and  non-breeding  adults  is  the  same,  which  implies  the  constraint  0 21  =  4>  12  +  022  (Model 
1.3).  In  this  model,  both  the  number  of  parameters  and  the  rank  of  the  Jacobian  are  3, 
indicating  that  all  three  parameters  can  be  estimated. 

For  the  inter-birth  emigration  model  of  the  right  whale  without  constraints  (Model  2.1), 
the  rank  of  the  Jacobian  is  3,  while  the  number  of  parameters  is  5.  If  we  assume  that 


Table  3.1:  Parameter  constraints,  number  of  parameters,  and  the  rank  of  the  Jacobian 
matrix  for  the  15  temporary  emigration  models.  *  indicates  that  values  of  all  the  free 
parameters  are  estimable. 


Species 

Model 

Constraints  on 
survival  probabilities 

Parameters  set 
to  known  values 

Number  of 
parameters 

Rank  of 
Jacobian 

Albatross 

1.1 

none 

none 

4 

3 

(Fig.  3.1b) 

1.2 

none 

Pi 

3 

2 

1.3* 

021  =  012  +  022 

none 

3 

3 

Right  Whale 

2.1 

none 

none 

5 

3 

(Fig.  3.1c) 

2.2 

none 

Pi 

4 

2 

2.3* 

021  =  032  =  013  +  033 

none 

3 

3 

Grey  Seal 

3.1 

none 

none 

6 

3 

(Fig.  3.2b) 

3.2* 

010  =  021  =  032 

=  (033  +  043)  =  044 

none 

3 

3 

3.3* 

010  =  021  =  032 

=  (033  +  043) 

Pi 

3 

3 

3.4* 

021  =  032 

=  (033  +  043)  =  044 

Pi 

3 

3 

Albatross 

4.1 

none 

none 

7 

3 

(Fig.  3.2c) 

4.2* 

010  =  021  =  032  =  043 

=  (044  +  054)  =  055 

none 

3 

3 

4.3* 

010  =  021  =  032  =  043 

=  (044  +  054) 

P5 

3 

3 

4.4* 

021  =  032  =  043 

=  (044  +  054)  =  055 

P5 

3 

3 

4.5* 

Logistic  Model1 

Pb 

3 

3 

1:  See  Equations  (3.10)-(3.17) 
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survival  probability  is  unaffected  by  reproduction  (<£21  =  </>32  =  </>i3  +  <^33;  Model  2.3),  both 
the  rank  of  the  Jacobian  and  the  number  of  parameters  become  3,  indicating  that  all  the 
three  remaining  parameters  (</>  13,  ^33,  and  p\ )  can  be  estimated. 

For  the  immature  emigration  model  of  the  grey  seal  without  constraints  (Model  3.1),  the 
rank  of  the  Jacobian  is  3,  while  the  number  of  parameters  is  6.  While  keeping  the  rank  of 
the  Jacobian  at  3,  we  can  reduce  the  number  of  parameters  to  3  in  several  ways  depending 
on  assumptions  on  how  mortality  varies  among  the  emigrated  stages.  Here  we  show  some 
specific  examples.  Model  3.2  assumes  that  survival  probabilities  of  all  stages  are  the  same. 
Under  this  assumption  the  three  parameters  ^>33,  <j> 43,  and  P4  can  be  estimated  separately, 
which  also  provides  estimates  of  the  survival  of  other  stages  (4> 21,  4> 32,  and  <^44).  The  next 
two  models  use  an  independent  estimate  of  the  capture  probability  of  stage  4  (jq).  This 
permits  estimation  of  an  extra  survival  probability.  For  example,  Model  3.3  assumes  that 
all  immature  stages  (0,  1,  2,  and  3)  have  the  same  survival  probability,  which  may  differ 
from  that  of  adults  (stage  4).  Model  3.4  assumes  that  the  survival  of  newborn  individuals 
(stage  0)  may  differ  from  that  of  all  other  stages. 

The  capture  probability  P4  required  for  Models  3.3  and  3.4  could  be  obtained  using 
Pollock’s  robust  design  (Pollock,  1982)  if  multiple  samples  of  the  mature  individuals  are 
collected  within  each  primary  sampling  occasion,  or  by  using  a  Cormack-Jolly-Seber  type 
mark-recapture  method  (e.g.  Burnham  et  al.,  1987;  Lebreton  et  ah,  1992)  if  capture  histories 
of  the  mature  individuals  over  multiple  sampling  occasions  are  available.  The  latter  method 
is  used  by  Clobert  et  ah  (1994)  to  estimate  survival  probability  and  age-specific  breeding 
probability.  This  method  also  provides  a  separate  survival  probability  estimate  for  mature 
stage  (<^44).  This  survival  probability  could  also  be  provided  in  these  models,  reducing  both 
the  number  of  parameters  and  the  rank  of  the  Jacobian  by  one. 

The  situation  for  the  immature  temporary  emigration  models  for  Wandering  Albatross 
is  the  same  as  the  grey  seal  models.  After  reducing  the  number  of  parameters  by  assuming 
that  survival  probability  of  all  stages  are  the  same  (<f> 21  =  4> 32  =  </>43  =  <l> 44  +  4>m  =  ^55; 
Model  4.2),  we  can  estimate  all  remaining  parameters  (<^>44,  $54,  ps)-  Alternatively,  different 
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distributions  of  survival  probability  during  temporary  emigration  can  be  used.  Model  4.3 
assumes  that  survival  probability  is  the  same  for  the  first  4  stages  (02i  =  032  =  043  = 
044  +  054)  and  a  separately  estimated  capture  probability  (ps)  is  provided.  Then  we  can 
estimate  the  remaining  three  parameters  (044,  054>  and  0 55)-  Similarly,  model  4.4  assumes 
that  survival  probability  is  the  same  for  the  last  4  stages  (0 32  =  043  =  044  +  054  =  05s)  and 
P5  is  provided.  Then  we  can  estimate  0 21,  044)  and  054,  separately.  This  result  suggests 
that  the  length  of  temporary  emigration  in  our  immature  emigration  model  does  not  change 
the  number  and  types  of  parameters  that  can  be  estimated. 

Another  approach  to  constraining  parameters  in  the  immature  temporary  emigration 
models  is  to  write  the  survival  probability  as  a  parametric  function  of  the  stage  and  estimate 
the  parameters  in  that  function.  For  example,  with  the  stage  structure  of  Wandering 
Albatross,  we  might  model  the  survival  probability  as  a  logistic  function  of  age  with  the  last 
immature  stage  and  mature  stage  having  the  same  survival  probability.  In  many  organisms, 
survival  probability  increases  with  age;  therefore,  this  model  may  be  realistic  in  many  cases. 
Let  Si  be  survival  probability  of  individuals  in  stage  i,  then 


exp  (7  +  5i ) 

1  +  exp  (7  +  5i) 


for  i  —  0, ...,  4, 


(3.10) 


where  7  and  5  are  intercept  and  slope  parameters.  From  the  survival  probability,  transition 
probabilities  can  be  calculated  as 


001 

=  so 

(3.11) 

021 

=  Si 

(3.12) 

• 

(3.13) 

043 

=  S3 

(3.14) 

044 

=  (1  0)s4 

(3.15) 

054 

=  0S4 

(3.16) 

055 

=  s5, 

(3.17) 
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where  ip  is  the  probability  of  transition  from  stage  4  to  stage  5  conditional  on  survival.  If 
p5  is  known,  the  rank  of  the  Jacobian  is  3,  permitting  estimation  of  7,  <5,  and  tp  separately. 

3.5  Bias  caused  by  temporary  emigration 

To  examine  the  bias  of  the  estimates,  we  applied  Models  1.3,  2.3,  3.4,  4.2,  and  4.5,  shown  in 
Table  1,  to  simulated  data.  Each  simulated  dataset  contained  75  individuals  marked  at  the 
first  sampling  occasion.  For  each  individual,  the  stage  at  the  next  sampling  occasion  was 
selected  at  random  from  the  possible  states  (including  death)  with  probabilities  given  by 
the  column  of  the  transition  matrix  corresponding  to  the  current  stage.  An  individual  in  a 
capturable  stage  was  captured  randomly  with  a  probability  equal  to  the  capture  probability 
of  that  stage.  For  Models  1.3  and  2.3,  the  data  sets  consisted  of  9  resampling  occasions;  for 
Models  3.4,  4.2,  and  4.5,  sequences  of  19  resampling  occasions  were  generated. 

From  the  data  we  estimated  the  parameters  in  the  corresponding  model.  This  process 
was  repeated  to  obtain  1000  sets  of  parameter  estimates.  From  the  same  data,  we  also 
estimated  survival  probability  under  the  erroneous  assumption  that  there  was  no  temporary 
emigration,  using  the  Cormack-Jolly-Seber  method  (Cormack,  1964;  Jolly,  1965;  Seber, 
1965)  with  constant  survival  and  capture  probabilities. 

Table  3.2  compares  the  values  of  parameters  that  were  used  to  generate  data  and  the 
means  of  1000  estimates  for  the  five  models.  Biases  in  both  capture  probability  and  transi¬ 
tion  probability  are  very  small.  Model  1.3  performs  most  poorly,  but  even  here,  the  survival 
probability  {<p 22  +  ^12)  appears  to  be  without  bias.  On  the  other  hand,  Table  3.3  compares 
actual  values  and  the  means  of  1000  estimates  under  the  false  assumption  that  there  was 
no  temporary  emigration.  Assuming  there  was  no  temporary  emigration  while  there  was 
in  fact  temporary  emigration  caused  overestimations  of  survival  probability  for  Models  1.3 
and  2.3  and  underestimations  for  most  stages  in  Models  3.4  and  4.2,  showing  that  these 
biases  are  a  serious  problem. 

In  the  parameter  estimation  processes  above,  we  were  not  always  able  to  estimate  pa- 
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rameters.  Under  Models  1.3  and  2.3  with  75  marked  individuals,  about  14%  and  9%, 
respectively,  of  the  numerical  optimization  failed  to  converge.  These  appear  to  be  primar- 


Table  3.2:  Actual  and  mean  of  estimated  parameters  for  Models  1.3,  2.3,  3.4,  4.2,  and  4.5. 
Values  in  parentheses  are  standard  deviations.  Results  based  on  simulating  1000  replicate 
datasets. 


Model 

Parameter 

Actual  Value 

Mean  Estimate  (SD) 

Model  1.3 

Pi 

0.80 

0.79  (0.09) 

0  22 

0.40 

0.38  (0.08) 

012 

0.50 

0.52  (0.08) 

Model  2.3 

Pi 

0.90 

0.90(0.06) 

033 

0.30 

0.30  (0.06) 

013 

0.60 

0.60  (0.05) 

Model  3.4 

P4 

0.90 

known 

010 

0.50 

0.50  (0.10) 

033 

0.40 

0.40  (0.08) 

043 

0.50 

0.51  (0.08) 

Model  4.2 

P5 

0.80 

0.80  (0.02) 

044 

0.45 

0.45  (0.06) 

054 

0.50 

0.50  (0.06) 

Model  4.5 

P5 

0.80 

known 

7 

1.3863 

1.40  (0.36) 

6 

0.3895 

0.39  (0.13) 

0 

0.4737 

0.46  (0.07) 

Table  3.3:  Estimated  capture  and  survival  probabilities  wrongly  assuming  no  temporary 
emigration  (simple  Cormack-Jolly-Seber  model  with  constant  capture  and  survival  proba¬ 
bility).  Data  were  generated  assuming  Models  1.3,  2.3,  3.4,  and  4.2.  Values  in  parentheses 
are  standard  deviations.  Results  based  on  simulating  1000  replicate  datasets. 


Model 

Parameter 

Actual  Value 

Mean  Estimate  (SD) 

Model  1.3 

Capture  Probability  ( p ) 
Survival  Probability  (s) 

0.80  for  stage  1 
0.90 

0.20  (0.02) 

0.94  (0.02) 

Model  2.3 

Capture  Probability  (p) 
Survival  Probability  (s) 

0.90  for  stage  1 
0.90 

0.14  (0.01) 

0.98  (0.01) 

Model  3.4 

Capture  Probability  (p) 
Survival  Probability  (s) 

0.90  for  stage  4 
0.5  for  stage  0 
0.9  for  stage  >  0 

0.52  (0.05) 

0.83  (0.03) 

Model  4.2 

Capture  Probability  (p) 
Survival  Probability  (s) 

0.80  for  stage  5 
0.95 

0.63  (<0.01) 

0.72  (<0.01) 
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ily  caused  by  extrinsic  identifiability  problems  as  repeating  the  analysis  with  150  marked 
individuals  reduced  the  incidence  of  non-convergence  to  5%  and  1%,  respectively. 

3.6  An  example:  the  right  whale 

As  an  example,  we  apply  our  method  to  data  on  the  North  Atlantic  right  whale  ( Eubalaena 
glacialis).  This  is  one  of  the  most  endangered  mammal  populations  in  the  world,  currently 
consisting  of  fewer  than  300  individuals.  There  is  evidence  to  suggest  that  this  number  is 
slowly  declining  because  of  declined  survival  probability  (Caswell  et  ah,  1999).  The  right 
whale  is  found  along  the  east  coast  of  the  United  States.  In  summer,  they  feed  in  the 
northwest  Atlantic,  including  Massachusetts  Bay,  Great  South  Channel,  Bay  of  Fundy,  and 
Brown’s  Bank.  In  winter,  some  females  migrate  to  the  coast  of  Florida  and  Georgia  for 
calving.  Their  inter-birth  interval  is  3-5  years. 

The  capture  history  data  we  use  here  are  based  on  photographs  taken  by  members  of 
the  North  Atlantic  right  whale  consortium  during  annual  surveys  along  the  east  coast  of  the 
United  States.  Because  right  whales  have  unique  markings  called  callosity  patterns  on  their 
heads,  individuals  can  be  identified  from  photographs.  Based  on  these  photographs,  the 
New  England  Aquarium  has  been  accumulating  capture  history  data  of  individuals  since 
1980.  Here,  we  use  only  sightings  of  mothers  (i.e.  females  with  their  calf).  We  consider 
individuals  to  have  been  marked  on  the  occasion  of  their  first  identification  as  mothers,  and 
recaptured  when  they  were  resighted  with  a  calf  during  a  subsequent  year.  If  the  data  had 
actually  been  collected  this  way,  females  would  exhibit  an  inter-birth  temporary  emigration 
with  an  interval  of  at  least  3  years.  This  corresponds  to  the  data  structure  in  the  ongoing 
study  of  southern  right  whales  ( Eubalaena  australis )  in  their  breeding  ground  off  Peninsula 
Valdes,  Argentina  (Payne  et  ah,  1990). 

We  use  the  stage  structure  shown  in  Figure  3.1c.  This  stage  structure  is  consistent  with 
the  fact  that  right  whales  do  not  give  birth  for  at  least  two  years  after  successful  repro¬ 
duction.  Because  previous  analyses  have  detected  time- variation  in  survival  and  transition 
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probabilities  (Caswell  et  al.,  1999),  we  divided  the  data  into  two  periods,  1980-1988  (t  =  1) 
and  1989-1997  (t  =  2).  Transition  and  survival  probabilities  are  assumed  constant  within 
each  period.  We  further  assumed  that  the  survival  probabilities  are  the  same  for  all  stages 
(Model  2.3).  The  Jacobian  indicates  that  this  model,  given  the  available  data,  is  identifi¬ 
able.  Because  individuals  in  stage  3  have  three  possible  fates  (remaining  in  stage  3,  moving 
to  stage  1,  and  death)  the  probability  of  which  must  sum  to  1,  the  transition  probabilities 
(013  and  0 33)  in  each  period  (t  =  1,2)  were  modeled  with  polychotomous  logistic  equations: 


At)  _  exp  {at) 

13  1  +  exp  (at)  +  exp  (f 3t ) 

At)  _  exp  ((3t) 

33  1  +  exp  (at)  +  exp  ($)  ’ 


(3.18) 

(3.19) 


where  at  and  fit  are  the  parameters  to  be  estimated.  Because  the  model  assumes  that  the 
survival  probability  of  all  stages  are  the  same,  0^  =  032  =  0i3  +  0 33  •  We  expressed  the 
likelihood  associated  with  the  data  using  the  method  in  Chapter  2,  and  the  negative  log 
of  the  likelihood  was  numerically  minimized  using  Matlab  function  fminu  to  estimate  the 
six  parameters  (pi,  P2,  «i,  a2,  Pi,  and  @2)- 

To  compare  the  estimated  parameters  between  the  two  periods  in  terms  of  demographic 
context,  we  calculated  survival  probability  as 


S(t)  =  013  +  033- 


(3.20) 


We  also  calculated  the  number  of  reproductive  events  based  on  Markov  chain  theory: 


Rw  =  F  (i  -  A®) 


(3.21) 
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where 


A«  = 


F  = 


'  8?  0  0  ' 

8?  0  0 

^  0  82  83  ) 

( 1  0  0  ^ 

0  0  0 
0  0  0  J 


(3.22) 


(3.23) 


and  I  is  the  identity  matrix.  The  value  r*n  gives  the  expected  number  of  future  reproductive 
events  from  a  female  that  has  just  given  a  birth.  More  details  of  this  calculation  can  be 
found  in  Chapter  2. 

Figure  3.3  compares  survival  probability  and  expected  number  of  reproductive  events 
between  the  two  periods.  Error  bars  in  the  figures  are  point-wise  95%  confidence  intervals 
based  on  a  parametric  bootstrap  sampling  procedure.  We  drew  bootstrap  samples  of  the 
four  parameters  5^,  S^2),  / and  /u2)  from  a  multivariate  normal  distribution,  with 
mean  vector  equal  to  the  maximum  likelihood  estimates  and  covariance  matrix  calculated 
from  the  inverse  of  the  last  Hessian  matrix  obtained  during  model  fitting.  A  thousand 
bootstrap  samples  were  generated  and,  for  each  bootstrap  sample  of  the  parameters,  the 
survival  probabilities  and  expected  number  of  reproductive  events  were  calculated.  The  95 
%  confidence  intervals  for  these  parameters  were  defined  by  the  2.5  and  97.5  percentiles  of 
the  simulated  distributions. 

As  expected,  Figure  3.3  shows  confidence  intervals  for  survival  and  reproductive  events 
are  wide  because  of  the  small  sample  size.  However,  the  results  suggest  that  both  survival 
and  the  expected  number  of  reproductive  events  have  declined.  This  result  is  consistent 
with  the  previous  finding  in  Fujiwara  &  Caswell  (2001)  although  the  declines  revealed  by 
the  analysis  of  the  full  data  set  are  greater  than  those  found  here. 
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Figure  3.3:  Analyses  of  North  Atlantic  right  whale  data,  (a)  Survival  probability  of  females 
that  have  previously  given  birth  at  least  once  during  periods  between  1980  and  1988  and 
between  1989  and  1997.  (b)  Expected  number  of  future  reproductions  during  the  lifetime 
of  females  that  have  given  birth  at  least  once. 

3.7  Discussion 

Temporary  emigration  causes  an  extreme  case  of  heterogeneity  in  capture  probability.  Tem¬ 
porarily  emigrated  individuals  have  zero  capture  probability,  while  the  rest  of  individuals 
have  a  non-zero  capture  probability.  Ignoring  temporary  emigration  when  it  exists  leads  to 
biased  estimates  of  survival  probability.  By  explicitly  modeling  the  emigration  process,  the 
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method  presented  in  this  paper  eliminates  this  bias. 

The  temporary  emigration  problem  is  often  encountered  in  sampling  large  animals  such 
as  seabirds,  sea  turtles,  and  marine  mammals.  When  the  same  problem  is  encountered  with 
smaller  organisms  such  as  fish  or  insects,  researchers  may  be  able  to  eliminate  the  temporary 
emigration  problem  by  experimental  manipulation.  For  example,  it  has  been  suggested  that 
when  organisms  can  be  marked  and  released  into  unobservable  stages  in  the  inter-birth 
temporary  emigration  stage  structures  (e.g.,  Models  3.1  and  4.1),  all  the  parameters  may 
be  estimated  separately  (G.  White,  personal  communication ).  In  such  cases,  the  methods 
outlined  here  to  determine  the  estimability  of  parameters  are  still  useful. 

Our  results  are  complementary  to  the  approach  of  Kendall  et  al.  (1997)  and  Kendall 
(1999),  who  use  the  Pollock’s  robust  design  method  (Pollock,  1982)  to  deal  with  tempo¬ 
rary  emigration.  In  their  method,  a  closed  population  model  is  used  to  estimate  capture 
probabilities,  which  are  incorporated  into  an  open  population  model  to  estimate  survival 
probabilities.  This  method  has  an  advantage  over  our  method  in  that  it  permits  estimation 
of  survival  probabilities  that  vary  freely  from  one  sampling  occasion  to  the  next.  The  ro¬ 
bust  design  method,  however,  requires  multiple  secondary  samplings  within  each  primary 
sampling  period  in  order  to  obtain  information  on  the  capture  probability.  The  method 
presented  here  does  not  require  secondary  sampling,  although  it  can  be  combined  with  the 
robust  design  method  when  such  samples  are  available. 

In  this  paper,  we  focused  on  two  stage  structures,  representing  two  different  emigration 
processes,  and  considered  only  selected  constraints  on  the  parameters.  We  emphasize  that 
the  method  can  be  applied  to  any  biologically  interesting  age-  or  stage-classified  life  cycle, 
with  any  pattern  of  observed  and  unobserved  stages.  Each  such  structure  may  admit  several 
choices  for  constraining  parameters  to  make  estimation  possible.  Choosing  biologically 
relevant  constraints  is  part  of  the  data  analysis  process.  We  are  certain  that  many  useful 
stage  structures  and  parameter  constraints  exist  and  are  yet  to  be  discovered. 
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Chapter  4 


Estimating  dispersal  kernels  from 
individual  mark-recapture  data 


4.1  Introduction 

A  dispersal  kernel  is  a  probability  density  function  for  the  location  of  an  individual  at  time 
t  +  1  (t  =  0, 1, 2, . . .)  conditional  on  the  location  of  the  individual  at  time  t.  This  kernel 
is  often  incorporated  in  an  integrodifference  equation  to  model  movements  of  individual 
organisms  (e.g.  Kot  et  al.,  1996;  Van  Kirk  &  Levis,  1996;  Van  Kirk  &;  Lewis,  1999;  Lewis  & 
Pacala,  2000;  Neubert  &  Caswell,  2000;  Neubert  et  al.,  2000;  Clark  et  al.,  2001b).  When  the 
kernel  is  the  probability  density  function  of  the  Normal  distribution  (Gaussian  distribution), 
the  corresponding  integrodifference  equation  model  is  a  discrete-time  equivalent  of  a  simple 
diffusion  process.  The  dispersal  kernel  can  take  other  shapes  to  represent  different  hypothe¬ 
ses  about  movements  of  individuals  (Neubert  et  al.,  1995).  For  example,  organisms  may 
tend  to  either  stay  where  they  are  or  to  travel  long  distances  when  they  decide  to  travel. 
Such  a  movement  pattern  can  be  modeled  with  leptokurtic  probability  density  functions 
such  as  the  Laplace  and  Cauchy  density  functions.  In  fact,  leptokurtosis  is  a  ubiquitous 
feature  of  dispersal  data  in  ecology  (Okubo,  1980;  Kot  et  al.,  1996). 

Recent  studies  in  theoretical  ecology  have  shown  that  different  movement  patterns  (i.e. 
different  shapes  of  the  kernels)  can  produce  qualitatively  different  population  dynamics. 
For  example,  Kot  et  al.  (1996)  showed  that  invading  organisms  with  frequent  long-distance 
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dispersals  (i.e.  “fat-tailed”  dispersal  kernels)  can  have  an  accelerating  speed  of  invasion. 
Clark  (1998)  showed  that  a  dispersal  pattern  modeled  by  the  fat-tailed  kernel  can  explain 
the  extremely  fast  spread  of  trees  during  the  Pleistocene.  In  other  studies,  Kot  (1989)  and 
Neubert  et  al.  (1995)  showed  that  dispersal-driven  instability  may  produce  stable  spatial 
patterns  of  population  density  in  predator-prey  systems  under  some  patterns  of  predator 
and  prey  dispersal.  Van  Kirk  &  Levis  (1996)  also  showed  that  the  scale  of  a  dispersal  kernel 
relative  to  the  size  of  a  fragmented  habitat  can  determine  persistence  of  a  population.  These 
theoretical  developments  are  prompting  ecologists  to  measure  dispersal  kernels  in  the  field. 

Here,  we  present  a  method  for  estimating  dispersal  kernels  from  individual  mark-recapture 
data.  The  mark-recapture  data  consist  of  information  on  whether  or  not  each  uniquely 
marked  individual  was  captured  at  each  sampling  occasion  along  with  its  captured  loca¬ 
tion.  From  the  data,  we  derive  a  likelihood  function,  and  a  maximum  likelihood  method  is 
used  to  estimate  parameters  associated  with  the  dispersal  kernels. 

When  an  individual  is  not  recaptured,  we  do  not  know  its  location.  One  possibility  is 
that  it  was  outside  of  the  sampling  area;  in  this  case,  there  was  zero  probability  of  it  being 
captured.  The  other  possibility  is  that  the  individual  was  actually  within  the  sampling  area 
but  was  not  captured  by  chance.  To  account  for  these  two  possibilities,  we  explicitly  model 
the  sampling  process  using  a  function  for  the  capture  probability  at  a  given  location;  we  call 
this  function  the  capture  probability  function.  The  capture  probability  function  is  0  where 
sampling  was  not  conducted,  and  bounded  between  0  and  1  where  sampling  was  conducted. 
The  capture  probability  function  and  dispersal  kernel  are  used  in  integrodifference  equations 
to  express  the  likelihood. 

In  Section  4.2  and  4.3,  we  describe  examples  of  dispersal  kernels  and  capture  probability 
functions,  respectively.  In  Section  4.4,  an  algorithm  for  formulating  the  likelihood  function 
is  provided.  In  Section  4.5,  we  show  examples  of  joint  probability  densities  of  locations  of 
individuals  at  time  t  and  their  capture  history  sequences  up  to  time  t  —  1.  These  probability 
densities  are  used  to  demonstrate  how  a  sampling  process  changes  the  probability  densities 
of  locations  of  individuals.  In  Section  4.6,  we  evaluate  the  performance  of  the  estimators 
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using  simulated  data  sets.  In  Section  4.7,  we  compare  some  results  in  Section  4.6  and 
equivalent  results  from  the  same  data  but  falsely  assuming  that  a  sampling  domain  was 
large  enough  to  include  all  individuals  at  all  sampling  occasions  and  falsely  assuming  a 
wrong  dispersal  kernel  shape.  In  Section  4.8,  we  examine  the  performance  of  the  method 
in  selecting  the  true  dispersal  kernel  shape  based  on  Akaike  Information  Criteria  (AIC). 
Up  through  Section  4.8,  we  assume  individuals  do  not  die,  or  that  deaths  occur  on  a  much 
longer  time-scale  than  spatial  movements.  In  Section  4.9,  we  discuss  how  we  can  incorporate 
the  deaths  in  a  likelihood  formulation  when  data  for  a  longer  time-scale  are  also  available. 
Other  extensions  of  the  spatial  mark- recapture  method  are  also  discussed  in  the  last  section. 

4.2  Dispersal  Kernels 

Gaussian,  Laplace,  and  Cauchy  functions  are  used  in  our  analysis  as  dispersal  kernels;  these 
functions  are  denoted  by  g,  l,  and  c,  respectively.  These  three  kernels  differ  in  kurtosis  and 
fatness  of  their  tails,  which  are  important  characteristics  that  can  determine  qualitatively 
different  ecological  dynamics  ( e.g .  Kot  et  al.,  1996;  Neubert  et  al.,  1995;  Van  Kirk  &  Levis, 
1996).  The  Gaussian  distribution  is  (by  definition)  nonleptokurtic  and  has  exponentially 
bounded  tails.  The  Laplace  distribution  is  leptokurtic  and  has  exponentially  decaying 
tails.  The  Cauchy  distribution  is  leptokurtic  and  its  tails  are  not  exponentially  bounded. 
Although  we  focus  on  these  three  distributions,  the  approach  we  outline  in  this  paper  can 
be  extended  to  many  other  shapes  of  dispersal  kernels. 

Let  At  €  3?  be  a  random  variable  for  the  location  of  an  individual  at  time  t  =  0, 1, 2, . . . ,  T 
where  T  is  the  time  of  the  last  sampling.  Then  a  kernel  kd{xt+i\xt)  of  shape  d  (d  —  g,  l,  c ) 
is  the  probability  density  function  for  Xt+i  conditional  on  Xt  —  xt-  With  the  further  as¬ 
sumption  that  the  kernel  depends  only  on  the  relative  locations  of  an  individual  at  times  t 
and  t  +  1,  the  Gaussian  kernel  is 

kg(xt+i  -xt)  =  — 7^:  exp  +  2a2  ~  ’  (4,1) 
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where  a  >  0.  Similarly,  the  Laplace  kernel  is 

ki  (®t+i  -  xt)  =  exP  ,  (4.2) 

where  (3  >  0.  Finally,  the  Cauchy  kernel  is 

kc(xt+i  ~  xt)  =  ^7  1  +  ^ ^  |  ,  (4.3) 

where  7  >  0.  In  the  above  density  functions,  a,  (3  and  7  are  parameters  that  define  the 
scale  of  dispersal  distances  (displacements),  and  p,  is  a  parameter  that  defines  a  directional 
movement  of  the  individual  between  consecutive  sampling  occasions  (i.e.  a  shift  in  the 
location  of  the  center  of  a  kernel). 


4.3  Capture  Probability  Functions 

The  capture  probability  function  pt{xt)  gives  the  probability  of  capturing  an  individual 
conditional  on  its  location  at  time  t  is  xt.  This  function  depends  on  how  sampling  is 
conducted.  For  example,  the  sampling  may  be  done  uniformly  within  a  sampling  area 
[- d,d ’].  Under  this  sampling  protocol  (Mi),  the  capture  probability  function  is 


Vt{xt)  = 


ut  for  \xt\  <  d, 
0  for  \xt\  >  d. 


In  addition  to  the  first  protocol,  an  intense  sampling  area  [-r,r]  within  the  sampling 
domain  may  be  added.  Under  this  second  protocol  (M2),  the  capture  probability  function 

is 

1  for  \xt\  <  r, 

Pt{xt)  =  <  ut  for  r  <  |xf|  <  d,  (4-5) 


In  the  third  protocol  (M3),  sampling  is  only  done  outside  of  the  center  region: 


Pt(xt)  =  { 


0 


(  ut 


for 

for 


N  <  d, 

\xt\  >  d. 


(4.6) 


This  type  of  sampling  protocol  is  sometimes  used  to  obtain  information  on  the  shape  of  tail 
regions. 

In  the  all  three  sampling  models,  ut  is  the  parameter  to  be  estimated.  The  other  para¬ 
meters  ( d  and  r)  axe  sizes  of  a  sampling  domain  and  an  intense  sampling  region,  respectively. 
Because  they  are  determined  before  sampling,  they  are  known  constants. 

Although  we  focus  on  capture  probability  functions  (4.4)-(4.6)  for  simplicity,  other  func¬ 
tions  can  be  incorporated  in  the  method  outlined  in  this  paper.  In  Section  4.9,  we  discuss 
other  functions  that  may  be  useful. 


4.4  An  Algorithm  for  the  Likelihood 

Let  6  be  a  parameter  set  that  defines  the  dispersal  kernel  (the  scale  and  shift  parameters) 
and  the  capture  probability  function  ( ut ).  A  likelihood  h(0)  associated  with  the  zth  in¬ 
dividual  in  the  data  is  a  function  of  0  and  depends  on  the  capture  history  of  individual 
%. 

Let  Sx}1+j{xtlt  ztj+i,  •  •  •  ,xti+j)  be  the  joint  probability  density  function  for  the  locations 
X<  (t  =  t\, . . . , ii  +  j)  of  individual  i  (i  =  1, . . . ,  n)  and  for  its  capture  history  sequence  up 
to  time  t\  +  j.  We  assume  that  individual  i  was  first  captured  (i.e.  marked)  at  time  t\  and 
at  location  xtl .  The  likelihood  can  then  be  calculated  by  the  following  algorithm: 

1.  Obtain  the  probability  density  function  for  the  locations  Xtx  and  Xtl+\  of  individual 
i: 

StH-iOzt^ti+i)  =  kd{xt  1+1  -  Xtx)-  (4-7) 

2.  Obtain  the  joint  probability  density  function  of  the  locations  Xtl,  Xtl+i,  and  Xtl+2 


64 


of  individual  i.  If  individual  i  was  captured  at  location  Xtj+i  at  time  ti  +  1: 

^ti+2('*'h  i  ®ti+l>  ®ii+  2)  =  (®ti )  ®ii+l  )Pti+l (®ti  ’f"  l)^'d(^'ti+2  ®ti+l)'  (4-8) 

If  individual  f  was  not  captured  at  ti  +  1: 

st?+2(a:ii > ^ti+ij ®ti+2)  =  f  sp+i(-2:){l  —  Pti+i{z)'}^d(xti+2  ~  z)dz.  (4-9) 

J—o o 

3.  Obtain  the  joint  probability  density  function  of  the  locations  Xtx,  X^+i,  •  •  • ,  Xt1+t+ 1 
of  individual  i.  If  individual  i  was  captured  at  location  Xfj+t  at  time  t\  + 1: 

(4.10) 

If  individual  i  was  not  captured  at  t\  +  t : 

)  £fi+l)  •  •  •  )  3^1+t+l)  =  /  sii+t{z){^~Pti+t{z)}^d{^ti+t+l~  z)dz.  (4.11) 

J  —00 

4.  Repeat  step  3  until  the  time  of  the  final  sampling  occasion  T. 

5.  Calculate  the  probability  density  of  the  final  observation  and  reinterpret  it  as  a  like¬ 
lihood.  If  individual  i  was  captured  at  location  xt  at  time  T, 

£i{0)  =  s$\xT).  (4.12) 

If  individual  %  was  not  captured  at  time  T, 

£i(G)=  f  s$(z)(l-pT{z))dz.  (4.13) 

J  —  OO 

□ 

It  should  be  noted  that  equation  (4.11)  is  a  convolution  of  two  functions  and  can  there- 
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fore  be  calculated  quickly  using  a  fast  Fourier  transform  convolution  method. 

The  above  algorithm  is  repeated  for  all  marked  individuals.  Assuming  that  individ¬ 
ual  movements  are  independent  of  one  another,  the  likelihoods  associated  with  individual 
capture  history  sequences  can  be  multiplied  to  get  the  likelihood  of  the  data  including  all 
marked  individuals.  This  likelihood  can  be  numerically  maximized  to  estimate  the  parame¬ 
ters. 

4.5  Joint  Probability  Density  of  Location  and  Past  Capture 
History  Sequence 

To  calculate  the  likelihood  in  the  previous  section,  we  used  the  joint  probability  density 
function  s^(a:o, . . . ,  xt)  of  the  locations  of  individuals  from  time  0  to  t.  Here,  we  assume 
the  individual  was  marked  at  time  0.  In  this  section,  we  show  how  •  •  •  >  xt)  may 

change  with  time. 

As  examples,  we  use  two  marked  individuals  released  at  different  locations  at  t  =  0. 
We  simulated  movements  of  the  two  individuals  using  the  Gaussian  kernel  (a  =  1,  /x  =  0), 
and  their  capture  was  simulated  using  the  M*  sampling  protocol  with  probability  u  =  0.85 
within  d  =  [—1,1]  and  u  =  0  outside.  The  marginal  probability  density,  which  does  not 
depend  on  past  capture  histories,  can  also  be  calculated  by  convolving  the  kernels  at  each 
time  step.  We  have  plotted  the  joint  probability  density  function  s[z\xq,  . . . ,  ®t)  sliced 
at  actual  capture  history  sequence  up  to  time  t-1  (solid  lines)  and  a  marginal  probability 
density  function  of  location  (dotted  lines)  at  time  t  =  1, 2, 3  in  Figure  4.1. 

Figures  4.1a  and  4. Id  show  the  joint  probability  densities  of  the  individuals  at  f  =  1. 
They  are  just  Gaussian  probability  densities  and  the  same  as  the  marginal  probability 
densities  because  they  do  not  have  any  previous  sampling  occasions. 

Figures  4.1b  and  4.1e  show  the  joint  probability  densities  (solid  lines)  of  the  two  indi¬ 
viduals  at  t  =  2  when  the  individuals  were  not  captured  at  t  =  1.  Because  they  were  not 
captured  at  t  =  1,  the  probability  that  individuals  were  within  the  sampling  domain  at 


0.2 
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Figure  4.1:  Joint  probability  density  of  location  of  an  individual  and  its  past  capture  history 
sequence  (solid  line)  and  the  marginal  probability  density  of  location  of  the  individual 
(dashed  line).  We  simulated  movements  of  the  two  individuals  using  the  Gaussian  kernel 
(a  =  l,fi  =  0),  and  their  capture  was  simulated  using  the  Mi  sampling  protocol  with 
probability  u  =  0.85  within  d  =  [—1,1]  and  u  —  0  outside,  (a)  The  probability  densities 
at  t  =  1  when  an  individual  was  released  at  x  —  0  at  t  =  0.  (b)  The  probability  densities 
at  t  =  2  when  an  individual  was  released  at  x  =  0  at  t  =  0  but  not  recaptured  at  t  =  1. 
(c)  The  probability  densities  at  t  =  3  when  an  individual  was  released  at  x  =  0  at  t  —  0 
but  not  recaptured  at  t  =  1  and  t  =  2.  (d)-(f)  The  same  as  (a)-(c)  except  the  individual  is 
released  at  x  =  0.5  at  t  =  0.  Vertical  lines  indicate  the  edges  of  the  sampling  domain. 
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t  =  1  was  reduced.  Compared  with  the  marginal  probability  densities  (dotted  lines),  the 
shape  as  well  as  the  height  of  the  joint  probability  densities  are  changed  noticeably. 

Figures  4.1c  and  4. If  show  the  joint  probability  densities  of  the  two  individuals  at  t  =  3 
when  the  individuals  were  captured  at  neither  t  —  1  nor  t  =  2.  Because  the  individuals  were 
not  captured  in  two  previous  sampling  occasions,  the  degree  of  concavity  in  the  middle  of 
the  joint  probability  densities  are  magnified.  Comparison  between  the  joint  and  marginal 
probability  densities  show  how  the  shapes  of  the  two  probability  densities  depart  from  each 
other  because  of  the  sampling  process. 

For  the  infinite  integrals  in  equations  (4.11)  and  (4.13),  we  truncated  the  kernels  at 
x-min  =  —333  and  xmax  =  333,  renormalized  them  so  that  the  kernels  integrate  to  1,  and 
the  equations  were  integrated  over  a  domain  [—1000, 1000].  This  assumes  that  individuals 
can  only  travel  the  maximum  distance  of  333  in  one  time  step.  This  truncation  has  almost  no 
effect  on  the  calculation  because  the  probability  of  individual  traveling  beyond  the  distance 
333  according  to  the  Gaussian  kernel  (a  =  1,/j,  =  0)  is  infinitesimally  small. 

Comparisons  of  the  two  individuals  (compare  Fig  4.1b  and  4.1c  with  Fig  4.1e  and  4. If, 
respectively)  demonstrate  that  the  shape  of  the  joint  probability  density  function  depends 
on  where  an  individual  was  released.  For  example,  if  an  individual  was  released  near  the 
edge  of  the  sampling  domain,  it  is  more  likely  to  leave  the  sampling  domain  than  if  it  was 
released  in  the  middle  of  the  sampling  domain  before  the  next  sampling  occasion.  These 
differences  in  probability  influence  the  probability  distribution  of  location  of  an  individual 
at  the  subsequent  sampling  occasions.  These  results  show  that  although  the  movements  of 
the  individuals  are  modeled  with  a  memoryless  process,  the  calculation  of  s^(x)  depends 
on  the  past  capture  history  sequences.  This  makes  further  simplification  of  the  likelihood 
calculation  difficult. 


4.6  Bias  of  Parameter  Estimates 


We  demonstrate  performance  of  estimators  using  simulated  data  sets.  We  considered  7 
scenarios  that  differed  in  kernel  shapes,  time-dependence  of  capture  probability,  and  size  of 
sampling  domain;  these  scenarios  are  listed  in  Table  4.1.  Under  each  scenario,  we  generated 
data  and  estimated  parameters  from  the  data.  Each  data  set  consisted  of  100  individuals 
marked  and  released  at  the  origin  at  t  =  0  and  resampled  4  subsequent  sampling  occasions 
using  the  sampling  protocol  Mi .  This  was  repeated  50  times  to  obtain  means  and  standard 
errors  of  the  estimates. 

Scenarios  1,  2,  and  3  differ  in  their  kernel  shapes;  they  are  Gaussian,  Laplace,  and 
Cauchy,  respectively  (Table  4.1).  Here,  we  assume  that  movement  of  individuals  is  sym¬ 
metric  so  that  they  are  always  centered  at  the  location  of  the  previous  capture.  We  also 
assume  that  capture  probability  is  constant  over  time.  Under  these  scenarios,  there  are 
two  parameters  to  be  estimated:  the  scaling  parameter  for  the  kernel  (a,  (3,  or  7)  and  the 
capture  probability  within  the  sampling  domain  (u). 

Table  4.1:  Sampling  and  movement  scenarios  used  to  examine  parameter  estimation  per¬ 
formance. 


Scenario 

Kernel 

Capture  Probability 

Sampling  Domain 

Scenario  1 

Gaussian 
(/i  =  0,  a  —  1.0) 

Constant  (0.85) 

[-1.5, 1.5] 

Scenario  2 

Cauchy 

(m  =  0,  7  =  1-0) 

Constant  (0.85) 

[-1.5, 1.5] 

Scenario  3 

Laplace 

(fi  =  0,p=  1.0) 

Constant  (0.85) 

[-1.5, 1.5] 

Scenario  4 

Gaussian 
(fi  =  0.2,  a  =  1.0) 

Constant  (0.85) 

[-1.5, 1.5] 

Scenario  5 

Gaussian 
(fj,  =  0,  a  =  1.0) 

Time  Varying 
(0.8,  0.5,  0.7,  &  0.9  ) 

[-1.5, 1.5] 

Scenario  6 

Gaussian 
(n  =  0,  a  =  1.0) 

Constant  (0.75) 

[-1.5, 1.5] 

Scenario  7 

Gaussian 
(/j,  =  0,  a.  =  1.0) 

Constant  (0.85) 

[-1.0, 1.0] 
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Under  Scenario  4,  the  Gaussian  kernel  is  nonsymmetric  and  shifted  by  /x.  This  shift  is 
equivalent  to  an  advection  term  in  a  simple  diffusion-advection  model.  Under  this  scenario, 
there  are  three  parameters  (a,  u,  and  (i)  to  be  estimated. 

Under  Scenario  5,  the  simulated  capture  probability  is  varied  at  each  sampling  occasion. 
In  actual  mark-recapture  studies,  it  is  very  difficult  to  maintain  a  constant  level  of  sampling 
effort,  and  it  is  often  even  more  difficult  to  keep  the  same  sampling  conditions  over  multiple 
sampling  occasions.  Therefore,  the  changes  in  capture  probability  among  sampling  occasions 
are  often  expected.  Under  this  scenario,  we  have  five  parameters  (a,  u\,  112,  U3,  and  114)  to 
be  estimated. 

Scenario  6  and  7  are  the  same  as  Scenario  1,  but  use  lower  capture  probability  and 
reduced  sampling  domain  size,  respectively. 

For  all  kernel  shapes,  we  truncated  the  kernels  at  xmin  =  —300  and  Xmax  =  333, 
renormalized  them  so  that  the  kernels  integrate  to  1,  and  the  equations  were  integrated 
over  a  domain  [-1000, 1000]  to  calculate  the  infinite  integrals  in  equations  (4.11)  and  (4.13). 
This  process  slightly  changed  the  heights  of  the  kernels  and  makes  the  shape  tighter  in  the 
middle  for  a  kernel  with  fatter  tails  (e.g.  the  Cauchy  kernel).  We  expect  that  this  truncation 
of  kernels  causes  small  overestimations  of  the  scale  parameter  and  the  capture  probability. 

Table  4.2  shows  the  result  comparing  true  and  estimated  parameter  values.  Most  of  the 
estimates  agree  reasonably  well  with  the  true  values.  However,  with  such  a  small  sample  size 
(50  estimates  for  each  model),  some  estimates  deviate  from  the  true  values.  In  Scenario  2, 

/ 3  is  underestimated,  and  in  Scenarios  3  and  4,  7  and  a,  respectively,  are  overestimated.  In 
part,  these  small  deviations  can  be  attributed  to  errors  associated  with  the  truncation  and 
the  numerical  integration  of  kernels.  Furthermore,  likelihood  methods  do  not  necessarily 
give  unbiased  estimates.  In  the  next  section,  we  will  show  how  these  deviations  change  when 
we  make  the  incorrect  assumptions  that  the  sampling  was  done  in  a  larger  sampling  domain 
than  the  actual  size  and  that  the  shape  of  dispersal  displacements  were  approximated  by  a 
Gaussian  kernel  when  in  fact  come  from  a  leptokurtic  distribution. 


4.7  Bias  Caused  by  Assuming  Large  Sampling  Domain 


When  analyzing  mark-recapture  data,  one  might  be  tempted  to  assume  that  (1)  the  kernel 
is  approximately  Normal  and  that  (2)  the  sampling  domain  was  large  enough  to  include  all 
individuals  at  all  sampling  occasions.  If  the  first  assumption  were  true,  the  integrodiffer- 
ence  model  would  be  a  discrete-time  equivalent  of  a  simple  diffusion  model.  If  the  second 
assumption  were  true  (which  is  rarely  the  case)  and  sampling  was  done  uniformly  over  the 
sampling  domain,  then  the  algorithm  for  calculating  likelihood  could  be  simplified  because 
the  shape  of  the  joint  probability  function  in  Section  4.5  would  remain  the  same  regardless 
of  where  individuals  were  released  and  recaptured  previously.  When  both  assumptions  are 
true,  we  can  obtain  the  value  of  a  that  maximizes  the  likelihood,  which  is  expressed  using 


Table  4.2:  Comparison  between  true  and  estimated  parameter  values  (sample  size  50). 


Scenario 

Parameter 

True  Value 

Estimated  Value  (SE) 

Scenario  1 

a 

1.005  (0.010) 

u 

0.850  (0.005) 

Scenario  2 

(3 

1.00 

0.976  (0.015) 

u 

0.85 

0.833  (0.007) 

Scenario  3 

7 

1.00 

1.075  (0.020) 

u 

0.85 

0.862  (0.008) 

Scenario  4 

a 

1.00 

1.037  (0.012) 

V 

0.2 

0.193  (0.010) 

u 

0.85 

0.866  (0.006) 

Scenario  5 

a 

1.00 

0.993  (0.012) 

«i 

0.80 

0.794  (0.008) 

U2 

0.50 

0.505  (0.009) 

U3 

0.70 

0.704  (0.010) 

U4 

0.90 

0.886  (0.011) 

Scenario  6 

a 

0.995  (0.012) 

u 

0.743  (0.007) 

Scenario  7 

a 

1.000 

1.001  (0.022) 

u 

0.85 

0.843  (0.014) 
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the  algorithm  in  Section  4.4,  by 


a  — 


If  (Mi 

nh W’ 


(4.14) 


where  n  is  the  number  of  observed  dispersal  distances,  hi  is  the  length  of  the  ith  observed 
interval,  and  Vi  is  the  time  taken  to  travel  that  distance  (see  Appendix).  The  right  side  of 
equation  (4.14)  is  a  standard  deviation  of  the  length  of  the  observed  dispersal  intervals  in 
a  unit  time. 

To  demonstrate  how  falsely  assuming  a  large  sampling  domain  causes  bias  in  the  para¬ 
meter  estimation,  we  applied  the  equation  4.14  to  the  data  generated  under  Scenarios  1,  5, 
6  and  7  in  the  previous  section. 

Furthermore,  when  the  true  kernel  was  the  Laplace  or  the  Cauchy  (Scenario  2  and  3, 
respectively),  we  applied  the  algorithm  presented  in  Section  4.4  but  inappropriately  assumed 
a  large  sampling  domain  ( d  =  m)  to  the  data  generated  in  the  previous  section. 

In  the  results  of  these  calculations,  any  deviations  from  the  true  vales  in  addition  to 
the  deviations  found  in  the  previous  section  are  caused  by  the  wrong  assumption  on  the 
sampling  domain  size.  The  estimates  from  these  calculations  and  the  estimates  in  the 
previous  section,  which  set  the  sampling  domain  sizes  appropriately,  are  compared  in  Table 
4.3. 


Table  4.3:  Comparisons  of  scale  parameters  (a,  /3,  or  7)  when  appropriately  setting  the  size 
of  the  sampling  domain  and  when  wrongly  assuming  that  the  sampling  domain  was  large 
enough  to  include  all  individuals  at  all  sampling  occasions  (sample  size  50). 


Scenario 

True  Value 

Estimated  Value  (SE) 
using  the  kernel  method 

Estimated  Value  (SE) 
using  equation  (4.14) 

Scenario  1 

1.00 

1.005  (0.010) 

0.747  (0.004) 

Scenario  2 

1.00 

0.976  (0.015) 

0.566  (0.004) 

Scenario  3 

1.00 

1.075  (0.020) 

0.428  (0.004) 

Scenario  5 

1.00 

0.993  (0.012) 

0.728  (0.005) 

Scenario  6 

1.00 

0.995  (0.012) 

0.734  (0.005) 

Scenario  7 

1.00 

1.001  (0.022) 

0.561  (0.004) 
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Comparing  the  results  under  the  false  and  correct  assumptions,  it  is  clear  that  if  the 
sampling  domain  was  falsely  assumed  to  be  large,  the  scaling  parameters  (a,  /?,  and  7)  are 
underestimated  (Table  4.3).  This  results  from  the  fact  that  any  long-distance  movements 
were  not  observed  at  all  because  they  were  actually  outside  of  the  sampling  domain;  however, 
they  were  assumed  to  be  always  observable  under  the  false  assumption.  As  expected,  the 
bias  increases  as  the  size  of  the  actual  sampling  domain  is  reduced  from  d  =  1.5  in  Scenario 
1  to  d  =  1.0  in  Scenario  7  because  more  individuals  are  located  outside  of  the  sampling 
domain  when  it  is  smaller. 

Figures  4.2a,  4.3a,  and  4.4a  show  the  kernels  based  on  the  5 th  largest  and  the  5th 
smallest  of  the  50  parameter  values  estimated  in  the  previous  section  under  Scenario  1,  2, 
and  3,  respectively,  along  with  the  true  kernel  that  was  used  to  generate  the  data.  In  all 
three  cases,  the  true  kernel  is  within  the  upper  and  lower  estimates.  For  the  same  scenarios, 
Figures  4.2b,  4.3b,  and  4.4b  show  the  kernels  based  on  the  5 th  largest  and  the  5 th  smallest 
of  the  50  estimated  parameter  values  under  the  false  assumption  (large  sampling  domain 
size).  These  are  compared  with  the  true  kernel.  For  all  three  scenarios,  the  estimated 
kernels  are  tighter  than  the  true  kernel. 

These  deviations  from  the  true  kernels  are  expected  to  be  even  more  severe  if  the  true 
kernel  were  the  Laplace  or  the  Cauchy  but  was  assumed  to  be  Gaussian.  This  is  because 
more  individuals  are  lost  to  unsampled  tail  regions  when  true  kernels  are  leptokurtic,  so 
analyses  that  assume  they  are  observable  produce  bigger  errors.  To  demonstrate  this,  we 
applied  estimator  (4.14)  to  the  data  generated  under  Scenario  2  (Laplace  kernel)  and  under 
Scenario  3  (Cauchy  kernel).  In  doing  so,  we  are  falsely  assuming  that  the  kernel  was 
Gaussian  as  well  as  assuming  that  the  sampling  domain  was  large  enough  to  include  all 
individuals  at  all  sampling  occasions.  Again,  we  plot  the  kernels  based  on  the  5 th  largest 
and  the  5 th  smallest  parameter  values  along  with  the  true  kernel  (Fig.  4.3c  and  4.4c).  In 
both  cases,  under  the  wrong  assumptions,  more  individuals  were  falsely  expected  to  be 
nearer  to  the  center  than  the  true  distributions  (Fig.  4.3c,  and  4.4c). 

These  results  suggests  that  it  is  important  to  record  where  sampling  was  conducted  as 
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(a) 


Displacement 


(b) 


Displacement 


Figure  4.2:  Kernels  based  on  the  5th  largest  and  the  5th  smallest  of  the  50  estimated 
parameter  values  (solid  lines)  and  the  true  kernel  (Gaussian)  used  to  generate  the  data 
(solid  line  with  triangles),  (a)  Parameters  estimated  applying  Gaussian  kernel  and  the 
true  sampling  domain  size  (algorithm  in  Section  4.4).  (b)  Parameters  estimated  applying 
Gaussian  kernel  but  assuming  the  sampling  domain  size  was  large  enough  to  include  all 
individuals  at  all  sampling  occasions  (Equation  4.14).  See  the  text  in  Section  4.6  for  how 
data  were  generated. 
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Figure  4.3:  Kernels  based  on  the  5th  largest  and  the  5th  smallest  of  the  50  estimated  pa¬ 
rameter  values  (solid  lines)  and  the  true  kernel  (Laplace)  used  to  generate  the  data  (solid 
line  with  triangles),  (a)  Parameters  estimated  applying  Laplace  kernel  and  using  an  appro¬ 
priate  sampling  domain  size  (algorithm  in  Section  4.4).  (b)  Parameters  estimated  applying 
Laplace  kernel  but  using  a  larger  than  actual  sampling  domain  size  (d  =  m  in  the  algo¬ 
rithm  in  Section  4.4).  (c)  Parameters  estimated  applying  Gaussian  kernel  and  assuming  the 
sampling  domain  size  was  large  enough  to  include  all  individuals  at  all  sampling  occasions 
(4.14).  See  the  text  in  Section  4.6  for  how  data  were  generated. 


75 


Figure  4.4:  Kernels  based  on  the  5th  largest  and  the  5th  smallest  of  the  50  estimated  para¬ 
meter  values  (solid  lines)  and  the  true  kernel  (Cauchy)  used  to  generate  the  data  (solid  line 
with  triangles),  (a)  Parameters  estimated  applying  Cauchy  kernel  and  using  an  appropriate 
sampling  domain  size  (algorithm  in  Section  4.4).  (b)  Parameters  estimated  applying  Cauchy 
kernel  but  using  a  larger  than  actual  sampling  domain  size  (d=  m  in  the  algorithm  in  Sec¬ 
tion  4.4).  (c)  Parameters  estimated  applying  Gaussian  kernel  and  assuming  the  sampling 
domain  size  was  large  enough  to  include  all  individuals  at  all  sampling  occasions  (4.14).  See 
the  text  in  Section  4.6  for  how  data  were  generated. 
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a  part  of  the  mark-recapture  data  and  use  the  kernel  that  best  approximates  the  true  one. 


4.8  Dispersal  Kernel  Selection  Performance 

When  analyzing  data,  we  often  do  not  know  the  true  shape  of  the  kernel  in  advance.  Thus 
we  must  fit  several  different  kernels  and  try  to  select  the  one  closest  to  the  true  kernel,  a 
challenging  task.  Especially  when  different  shapes  of  kernels  result  in  different  population 
dynamics,  it  is  very  important  that  we  can  select  the  best-fit  kernel.  Here,  we  evaluate 
the  performance  of  Akaike’s  Information  Criteria  (AIC)  (Akaike,  1973)  as  a  tool  for  kernel 
selection. 

We  generated  a  data  set  using  each  of  the  three  dispersal  kernels,  fitted  the  three  disper¬ 
sal  kernels  to  the  data  set,  and  selected  the  best-fit  kernel  according  to  AIC.  Each  data  set 
consisted  of  100  individuals  released  at  the  origin  at  t  =  0  and  resampled  at  3  subsequent 
occasions.  This  was  repeated  20  times  to  obtain  the  probability  of  selecting  the  true  model. 
The  same  process  was  repeated  with  different  sampling  domains  ( d )  under  the  sampling 
protocol  Mi .  For  each  true  kernel,  we  set  the  scaling  parameter  (a,  /?,  or  7)  equal  to  \[2 
and  set  the  capture  probability  (u)  to  0.85  within  the  sampling  domain. 

Figures  4.5,  4.6,  and  4.7  show  changes  in  the  probability  of  selecting  the  true  kernel 
shape  with  the  size  of  sampling  domain  when  the  true  kernel  was  Gaussian,  Laplace,  and 
Cauchy,  respectively. 

Associated  error  bars  show  the  point-wise  estimates  of  95%  confidence  intervals.  In  each 
figure,  we  also  show  the  size  of  the  sampling  domain  that  contains  50%  of  individuals  at 
time  t  =  1  (left  line)  and  at  time  t  =  3  (right  line).  We  call  the  two  points  reference  points 
and  compare  the  probabilities  for  the  three  kernels  at  these  points. 

As  expected,  the  probability  of  selecting  the  true  kernel  shape  increases  as  the  sampling 
domain  increases  in  all  three  cases.  Comparisons  among  the  three  kernels  reveal  that  when 
the  true  kernel  was  Cauchy,  there  is  a  high  probability  of  selecting  the  true  kernel  even 
when  the  sampling  domain  is  small  (Figure  4.7).  When  the  true  kernel  was  Gaussian,  the 
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Size  of  Sampling  Domain 


Figure  4.5:  The  probability  of  selecting  the  true  kernel  shape  when  the  true  kernel  shape 
was  Gaussian  as  the  sampling  domain  is  changed.  Error  bars  indicate  estimated  point-wise 
95%  confidence  intervals.  Vertical  lines  on  the  left  and  right  indicate  the  size  of  sampling 
domain  where  50%  of  individuals  are  expected  to  be  within  the  sampling  domain  at  t  =  1 
and  t  —  3,  respectively. 
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Figure  4.6:  The  probability  of  selecting  the  true  kernel  shape  when  the  true  kernel  shape 
was  Laplace  as  the  sampling  domain  is  changed.  Error  bars  indicate  estimated  point- wise 
95%  confidence  intervals.  Vertical  lines  on  the  left  and  right  indicate  the  size  of  sampling 
domain  where  50%  of  individuals  are  expected  to  be  within  the  sampling  domain  at  t  =  1 
and  t  =  3,  respectively. 
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Figure  4.7:  The  probability  of  selecting  the  true  kernel  shape  when  the  true  kernel  shape 
was  Cauchy  as  the  sampling  domain  is  changed.  Error  bars  indicate  estimated  point- wise 
95%  confidence  intervals.  Vertical  lines  on  the  left  and  right  indicate  the  size  of  sampling 
domain  where  50%  of  individuals  are  expected  to  be  within  the  sampling  domain  at  t  =  1 
and  t  =  3,  respectively. 

probability  of  selecting  the  true  kernel  is  also  relatively  high  at  the  two  reference  points. 
The  probability  of  selecting  the  true  kernel  when  it  was  actually  Gaussian,  declines  rapidly 
to  the  left  of  the  smaller  reference  point.  We  speculate  that  this  rapid  decline  is  due  to  the 
nonleptokurtic  nature  of  the  distribution.  The  middle  part  of  the  kernel  tends  to  be  flat 
and  it  may  not  contain  much  information  about  the  shape  of  the  kernel.  When  the  true 
kernel  was  Laplace  function,  the  probability  of  selecting  the  true  kernel  is  the  lowest  (Figure 
4.7).  This  result  is  consistent  with  the  fact  that  the  shape  of  the  Laplace  distribution  is 
intermediate  between  Gaussian  and  Cauchy  distributions  in  terms  of  their  kurtosis. 

When  the  true  kernel  was  Cauchy  and  the  true  kernel  was  not  selected,  wrongly  selected 
kernels  tended  to  be  Laplace  kernel.  Similarly,  when  the  true  kernel  was  Gaussian,  wrongly 
selected  kernel  tended  to  be  Laplace.  However,  when  the  true  kernel  was  Laplace,  both 
Cauchy  and  Gaussian  kernels  were  selected  wrongly.  This  result  also  suggests  that  select¬ 
ing  between  Laplace  and  Cauchy  or  between  Laplace  and  Gaussian  is  more  difficult  than 
selecting  between  Cauchy  and  Gaussian.  If  such  model  selection  is  important,  one  should 
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consider  ways  to  improve  the  model  selection  performance.  Two  obvious  ways  to  improve 
the  performance  are  to  increase  sampling  domain  size  and  number  of  marked  individuals. 

Alternatively,  we  can  fix  the  size  of  sampling  domain  and  add  small  area  of  the  intensive 
sampling  within  the  domain  to  try  to  improve  chances  of  selecting  right  kernel.  We  examined 
how  the  probability  of  selecting  the  true  kernel  changes  as  the  size  of  intense  sampling  region 
is  changed  in  M2  sampling  protocol  (Fig.  4.8).  We  used  the  same  Gaussian  kernel  as  before 
(a  =  V 2  and  fx  =  0)  and  the  sampling  domain  d  =  0.75.  Without  the  intense  sampling 
domain,  there  is  only  0.45  probability  of  selecting  the  true  kernel.  However,  adding  a  small 
intense  sampling  region  r  =  0.1  in  Equation  (eq:intense)  increases  this  probability  to  0.8. 
This  suggests  that  intense  sampling  over  a  small  part  of  the  domain  provides  an  alternative 
way  of  improving  the  model  selection  performance  instead  of  simply  increasing  the  sampling 
domain  size  or  increasing  the  number  of  marked  individuals. 

Finally,  we  examined  the  performance  of  selecting  the  true  kernel  under  the  sampling 
protocol  M3,  in  which  only  outer  region  is  sampled.  This  protocol  may  be  used  when  trying 
to  collect  information  on  the  shape  of  the  tail  regions  of  kernels.  Under  this  protocol,  there 


Size  of  Intense  Sampling  Region 

Figure  4.8:  The  probability  of  selecting  the  true  kernel  shape  when  the  true  kernel  shape 
was  Gaussian  as  the  region  of  intense  sampling  was  increased  under  sampling  protocol  M2. 
Error  bars  indicate  estimated  point-wise  95%  confidence  intervals. 
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was  high  probability  of  selecting  the  true  kernel.  However,  when  the  sampling  region  was 
reduced  by  increasing  d  in  equation  4.6,  very  small  number  of  individuals  were  recaptured. 
This  was  especially  true  when  the  kernel  was  Gaussian.  As  a  result,  it  was  very  difficult 
for  likelihood  to  converge  in  the  numerical  optimization  even  when  applying  the  true  kernel 
shape.  This  suggests  that  although  the  important  part  of  the  kernel  in  terms  of  ecological 
dynamics  may  be  the  distribution  in  the  tail  regions,  it  is  also  necessary  to  sample  the 
center  region  to  estimate  kernels. 

4.9  Extension  of  the  Method 

We  have  considered  estimation  and  performance  of  selecting  kernel  shapes  under  different 
sampling  protocols;  however,  the  method  is  flexible  enough  to  be  used  in  other  circum¬ 
stances.  Here,  we  discuss  possible  extensions  of  the  method. 

4.9.1  Kernel  Shape 

We  used  three  kernel  shapes  (Gaussian,  Laplace,  and  Cauchy)  in  this  paper  because  of 
differences  in  their  kurtosis.  Other  kernel  shapes  can  be  incorporated  in  the  method,  but 
there  are  exceptions.  For  example,  a  rectangular  kernel  cannot  be  incorporated  in  the 
method.  The  parameters  in  the  rectangular  kernel  define  the  height  of  the  kernel  and 
the  parameters  of  the  capture  probability  function  also  define  the  height  of  the  function. 
As  a  result,  the  parameters  in  the  rectangular  kernel  and  capture  probability  function 
become  confounded.  It  should  also  be  noted  that  the  individual  mark-recapture  method 
permits  estimation  of  scale  parameters  (a,  /?,  and  7)  and  a  shift  parameter  (/x),  however, 
a  shape  parameter  such  as  the  one  in  the  standard  Student’s  t  distribution  appears  to  be 
unestimable.  The  estimated  parameters,  therefore,  are  conditional  on  the  shape  of  kernel, 
and  the  model  selection  approach  to  select  the  best-fit  kernel  as  presented  in  Section  4.8 
becomes  important. 

In  this  paper,  data  in  a  one-dimensional  space  were  considered.  However,  it  is  simple  to 
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extend  the  method  to  data  in  a  two-dimensional  space.  One  way  of  doing  so  is  to  collapse 
a  two-dimensional  into  a  one-dimensional  problem  using  a  polar  coordinate  if  the  kernels 
are  symmetric  about  their  centers  (see  Shigesada  &;  Kawasaki,  1997,  for  this  approach  in  a 
simple  diffusion  model).  Alternatively,  we  can  fit  two-dimensional  kernels  such  as  bivariate 
Gaussian  and  bivariate  Cauchy  distributions  directly  to  the  data. 

4.9.2  Mortality 

The  method  presented  here  assumes  that  there  were  very  few  deaths  of  marked  individu¬ 
als  within  the  entire  sampling  period  ( t  =  0, ...  ,T).  This  assumption  is  needed  because 
deaths  would  add  another  reason  for  not  capturing  individuals.  When  there  is  mortality, 
survival  probability  will  be  multiplied  to  the  dispersal  kernel.  The  product  is  the  kernel 
with  a  reduced  height.  A  kernel  with  a  similar  height  can  also  be  produced  by  letting 
more  individuals  to  travel  outside  of  sampled  region.  As  a  result,  kernel  parameters  become 
confounded  with  survival  probability.  The  same  problem  is  encountered  in  the  discrete 
multi-stage  mark-recapture  method,  in  which  individuals  temporarily  emigrate  out  of  sam¬ 
pled  region  (Kendall  et  al.,  1997;  Kendall,  1999,  Chapter  3).  In  the  discrete  case,  some 
transition  probability  parameters  are  also  confounded  (see  Chapter  3) . 

Under  certain  circumstances,  however,  parameters  associated  with  survival  probability 
can  be  estimated  along  with  the  dispersal  kernel  and  capture  probability.  One  such  example 
is  when  the  dispersal  kernel  has  finite  support  and  the  sampling  is  conducted  over  the  entire 
region  -within  which  individuals  move.  Then  there  is  no  temporary  emigration,  and  survival 
probability  can  be  estimated  in  addition  to  the  parameters  associated  with  the  kernel.  In 
this  case,  survival  probability  is  given  by  the  area  under  the  kernel. 

Alternatively  one  can  collect  data  on  two  different  time-scales  (one  on  the  scale  of 
movement,  and  the  other  on  the  scale  of  the  mortality  process)  to  estimate  parameters 
associated  with  both  kernel  and  survival  probability.  Movement  on  the  faster-time  scale 
can  be  modeled  assuming  no  mortality.  The  data  on  the  slower  time-scale  can  be  modeled 
with  dispersal  kernel,  survival  probability,  and  capture  probability.  A  single  likelihood 
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associated  with  both  sets  of  data  can  be  used  to  estimate  all  the  parameters.  This  method 
is  similar  to  Pollock’s  robust  design  method  (Pollock,  1982)  in  which  no  mortality  was 
assumed  in  a  short-time  scale  to  estimate  capture  probability. 


4.9.3  Sampling  Design 


In  this  paper,  we  considered  only  three  different  types  of  sampling  protocols.  However, 
there  are  many  other  possibilities.  In  Section  4.8,  we  demonstrated  that  even  with  a  small 
sampling  domain  relative  to  the  shape  of  the  kernel,  adding  a  small  intense  sampling  region 
within  the  sampling  domain  (sampling  protocol  M2)  can  improve  model  selection  perfor¬ 
mance.  For  the  protocol  we  considered  ,  we  added  the  intense  sampling  region  in  the 
middle  of  the  sampling  domain.  However,  the  intense  sampling  region  can  be  moved  to 
other  locations  within  the  sampling  domain.  It  is  also  possible  to  sample  at  many  isolated 
small  regions  intensively  or  also  possible  to  sample  only  toward  the  left  (or  right)  half  of 
the  space.  How  the  sampling  effort  should  be  allocated  is  an  important  question  to  ask 
before  actual  sampling  and  will  likely  vary  from  study  to  study.  To  answer  such  a  question, 
different  movement  and  sampling  scenarios  can  be  explored  using  simulated  data  sets  as  we 
demonstrated  in  this  paper. 

In  actual  sampling,  it  may  not  be  possible  to  maintain  a  constant  capture  probability 
over  the  whole  sampling  domain.  For  example,  vegetation  coverage  or  topography  of  the 
ground  may  influence  capturability  of  individuals.  However,  when  such  data  as  degree  of 
vegetation  coverage  or  index  of  topographic  feature  are  available,  they  can  be  incorporated 
into  the  capture  probability  function  as  covariates  using  a  function  such  as  the  logistic 


function 


Pt{xt)  = 


exp  (ip  +  u>c(xt)) 

1  +  exp  (ip  +  wc(ij)) 


(4.15) 


where  ip  and  1 0  are  the  intercept  and  slope  parameters  and  c(xt )  is  the  value  of  covariate  at 
location  x  at  time  t.  This  ability  to  model  the  capture  probability  function  is  one  strength 
of  the  individual  mark-recapture  method  presented  in  this  paper. 
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4.10  Discussion 


The  individual  mark-recapture  method  presented  here  permits  estimation  of  parameters 
associated  with  non-Gaussian  as  well  as  Gaussian  kernels  from  individual  mark-recapture 
data  that  are  sampled  in  a  finite  sampling  domain.  We  analyzed  the  performance  of  esti¬ 
mators  and  the  model  selection  method  using  simulated  data  sets.  Each  data  set  consisted 
of  only  100  marked  individuals,  but  both  the  parameter  estimators  and  the  model  selection 
worked  reasonably  well. 

The  method  presented  in  this  paper  can  be  considered  as  a  continuous-stage  equivalent  to 
a  multi-stage  (or  multi-type)  mark-recapture  method  (e.g.  Nichols  et  al.,  1992b,  also  Chapter 
2).  In  the  method  presented  in  this  paper,  “stage”  corresponds  to  the  spatial  location  of 
the  individual.  There  are  several  studies  that  attempt  to  quantify  movements  of  individuals 
using  a  discrete-stage  multi-stage  mark-recapture  method  (e.g.  Arnason,  1973;  Hilborn, 
1990;  Hestbeck  et  al.,  1991;  Brownie  et  al.,  1993).  However,  the  discrete-space  approach 
loses  the  detailed  information  on  where  individuals  were  located  within  the  discrete  patches. 
As  shown  in  Section  4.5,  it  is  very  important  information  when  estimating  dispersal  kernels 
accurately. 

Other  types  of  data  are  also  collected  to  quantify  movements  of  organisms.  For  example 
“seed  shadow”  data  contain  locations  of  seeds  released  from  a  single  plant.  From  these  data, 
non-Normal  dispersal  kernels  (e.g.  Werner,  1975;  Carey  &:  Watkinson,  1993)  or  invasion 
speed  of  the  population  (Clark  et  al.,  2001a)  can  be  estimated.  The  data  are  excellent  when 
organisms  have  both  dispersal  and  sessile  stages  and  the  source  of  seeds  can  be  identified. 
Another  type  of  data  is  mass  mark-recapture  data  in  which  markings  are  not  unique  to 
each  individual  (e.g.  Dobzhnsky  &  Wright,  1943;  Taylor,  1978;  Kareiva,  1983;  Turchin  & 
Thoeny,  1993).  These  data  contain  locations  of  captured  individuals,  but  do  not  contain 
information  on  which  individuals  were  captured.  This  type  of  data  is  excellent  when  large 
number  of  individuals  can  be  marked  and  sampling  process  can  be  controlled  to  be  constant 
over  space  and  time.  The  individual  mark-recapture  method  presented  in  this  paper  is  an 


addition  to  these  existing  ways  for  estimating  dispersal  kernel.  We  expect  this  method  to  be 
particularly  useful  to  estimate  kernels  associated  with  animals  that  can  potentially  travel  a 
long  distance  when  only  small  numbers  of  them  can  be  marked. 
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4.11  Appendix 

We  derive  equation  (4.14)  from  the  likelihood  function  in  Section  4.4.  Here,  we  are  assuming 
that  movement  is  modeled  with  the  Gaussian  kernel  without  any  shift  (p  —  0),  capture 
probability  is  constant  in  time  and  space,  and  sampling  is  done  over  the  entire  universe. 

First,  we  divide  the  capture  history  into  segments  when  the  individual  was  recaptured, 
and  note  that  the  probability  density  for  the  capture  history  sequence  is  a  product  of  the 
probability  densities  for  the  segments.  For  example,  suppose  that  an  individual  was  marked 
and  released  at  x\  at  time  1,  recaptured  the  first  time  at  x%  at  time  3,  but  never  recaptured 
again  when  the  sampling  ended  at  time  5.  This  capture  history  is  separated  into  two 
segments:  one  from  time  1  to  3  and  the  other  from  time  3  to  5.  The  probability  density  for 
the  second  segment  is  conditional  on  the  release  at  x%  at  time  3.  Consequently,  the  product 
of  the  probability  density  for  the  first  and  second  segments  gives  the  probability  density  of 
the  entire  capture  history  sequence. 

After  separating  the  capture  history  sequences  of  individuals  into  segments,  we  classify 
the  segments  into  two  forms.  In  one  form,  individuals  are  released  at  xt  at  time  t  and 
recaptured  at  xt+v  at  time  t  +  v  ( v  =  1,  2, . . .)  the  first  time  since  t,  and  in  the  other 
form,  individuals  are  released  at  xt  at  time  t  but  never  recaptured  again.  With  these  two 
forms,  we  can  construct  any  capture  history  sequences.  We  derive  general  expressions  for 
the  probability  density  of  these  two  types  of  the  segments. 

For  the  j th  segment  of  the  first  form,  we  write  down  the  probability  density  function 
for  the  location  Xt+V  =  y  of  an  individual  in  which  it  was  released  at  xt  at  time  t.  If  the 
individual  was  recaptured  immediately  after  the  release  (v  —  1), 


4+1  (y)  =  kg(y-xt)p 


(4.16) 
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where 


(4.17) 


kg{y  Xt)  -  exp  (y2a2  }  . 

Similarly,  if  the  individual  was  recaptured  the  first  time  at  t  +  2  (v  —  2), 


:+2(y)  =  [  kg{z  -  xt){l  -  p)kg(y  -  z)pdz 

J — OO 

roo 

=  (1  -p)p  kg(z  -  xt)kg(y  -  z)dz 

J— OO 


(4.18) 

(4.19) 


If  the  individual  was  recaptured  the  first  time  at  t  +  3  (v  =  3), 

, ..  roo  roo 

st+3(y)  =  /  /  fcg(z-xt)(l-p)fcff(fn-z)dz(l-p)fc5(j/-u;)pdu; 

J —OO  j  — OO 

/oo  roo 

/  A:5(2  -  xt)kg(w  -  z)dzkg(y  -  w)dw.  (4.20) 

-OO  J  —  OO 

By  inspection  of  the  above  equations,  we  obtain  the  general  probability  density  function 
associated  with  the  first  form  of  the  capture  history  segments  as 

4+v{y)  =  (1  -p)v~lpfv{y,xt),  (4.21) 

where  fv{y,xt )  is  the  probability  density  function  obtained  by  convolving  the  Gaussian 
kernel,  which  is  centered  at  xt,  with  v  —  1  Gaussian  kernels,  which  are  centered  at  the 
origin,  sequentially.  By  convolving  two  Gaussian  kernels,  we  get  the  Gaussian  kernel  with 
its  variance  being  the  sum  of  the  variances  of  the  two  convolved  Gaussian  kernels.  Therefore, 
equation  (4.21)  becomes 

=  t1  -  exp  [zwL]  ■  (422) 

Similarly  for  the  second  form  of  the  capture  history  sequence,  which  ends  without  ob¬ 
serving  the  individual,  we  write  down  the  probability  of  not  recapturing  the  individual  i 
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that  was  released  at  xt  at  time  t: 


s?  =  (1  ~p)  f  f\{z-xt)dz 

J  —oo 
r°° 

s?  =  (1  -p)2  /  h{z,xt)dz 
St  =  {l-pf  (  f3(z,xt)dz 

J —oo 


for 

T-t=  1, 

(4.23) 

for 

T-t  =  2, 

(4.24) 

for 

eo 

II 

1 

(4.25) 

where  T  is  the  time  of  the  last  sampling  occasion.  The  integrals  in  equations  (4.23)-(4.25) 
are  1  because  integrands  are  probability  density  functions.  Thus,  we  obtain  the  general 
equation  for  the  probability  as 

4?)  =  (l-pf-*.  (4.26) 


The  likelihood  in  Section  4.4  (Equations  4.12  and  4.13)  can  be  expressed  as  a  product  of 
equations  (4.22)  and  (4.26).  In  both  equations(4.22)  and  (4.26),  the  capture  probability  (p) 
can  be  factored  from  the  dispersal  kernels  because  the  capture  probability  is  independent 
of  location.  Therefore,  the  likelihood  can  be  expressed  by  a  product  of  two  terms  each  con¬ 
sisting  of  only  dispersal  kernel  or  only  capture  probability.  Consequently,  we  can  maximize 
the  terms  for  capture  probability  and  dispersal  kernel  the  likelihood  separately  to  find  the 
maximum  likelihood  estimates  of  the  parameters.  In  other  words,  capture  probability  does 
not  affect  the  parameters  associated  with  kernels  in  the  estimation.  Since  we  are  only  in¬ 
terested  in  estimating  the  parameter  associated  with  the  dispersal  kernel,  we  can  ignore  the 
capture  probability  part.  Furthermore,  equation  (4.26)  contains  only  the  capture  probabil¬ 
ities.  Thus  the  part  of  the  likelihood  that  we  are  interested  is  just  a  product  of  equations 
(4.22)  each  expressed  for  observed  dispersal  intervals.  We  denote  this  special  likelihood  as 
£sp.  After  collecting  terms  in  £sp,  we  obtain 
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where  n  is  the  number  of  observed  dispersal  intervals,  hj  is  the  length  of  the  jth  observed 
interval,  and  Vj  is  the  time  taken  to  travel  that  distance. 

Taking  the  derivative  of  equation  (4.27)  with  respect  to  the  parameter  a,  setting  the 
derivative  equal  to  zero,  and  solving  for  a  results  in  equation  (4.14). 
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Chapter  5 

Demography  of  the  endangered 
North  Atlantic  right  whale 
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Demography  of  the  endangered 
North  Atlantic  right  whale 

Masaml  Fujtwara  &  Hal  Caswell 

Biology  Department,  MS34,  Woods  Hole  Oceanographic  Institution,  Woods  Hole, 
Massachusetts,  USA 

Northern  right  whales  (Eubalaena  gladalis)  were  formerly  abun¬ 
dant  in  the  northwestern  Atlantic,  but  by  1900  they  had  been 
hunted  to  near  extinction.  After  the  end  of  commercial  whaling 
the  population  was  thought  to  be  recovering  slowly;  however, 
evidence  indicates  that  it  has  been  declining  since  about  1990 
(ref.  1).  There  are  now  fewer  than  300  individuals,  and  the  species 
may  already  be  functionally  extinct2,3  owing  to  demographic 
stochasticity  or  the  difficulty  of  females  locating  mates  in  the 
vast  Atlantic  Ocean  (Allee  effect4).  Using  a  data  set  containing  over 
10,000  sightings  of  photographically  identified  individuals  we 
estimated  trends  in  right  whale  demographic  parameters  since 
1980.  Here  we  construct,  using  these  estimates,  matrix  population 
models  allowing  us  to  analyse  the  causes  of  right  whale  imperil- 
ment.  Mortality  has  increased,  especially  among  mother  whales, 
causing  declines  in  population  growth  rate,  life  expectancy  and 
the  mean  lifetime  number  of  reproductive  events  between  the 
period  1980-1995.  Increased  mortality  of  mother  whales  can 
explain  the  declining  population  size,  suggesting  that  the  popula¬ 
tion  is  not  doomed  to  extinction  as  a  result  of  the  Allee  effect.  An 
analysis  of  extinction  time  shows  that  demographic  stochasticity 
has  only  a  small  effect,  but  preventing  the  deaths  of  only  two 
female  right  whales  per  year  would  increase  the  population 
growth  rate  to  replacement  level. 

Conservation  biology  uses  population  models  to  assess  popula¬ 
tion  performance,  to  diagnose  the  causes  of  poor  performance,  to 
prescribe  management  interventions,  and  to  make  prognoses  of 
population  viability  (see  chapter  18  of  ref.  5).  We  have  developed  a 
stage- structured  matrix  population  model5-7  that  addresses  all  four 
of  these  tasks.  This  matrix  model  uses  the  life  cycle  shown  in  Fig.  1  (a 
similar  model  has  been  successfully  applied  to  a  killer  whale 


Figure  1  Life  cycle  graphs  of  female  (a)  and  male  (b)  right  whales.  Numbers  represent 
different  stages:  1 ,  calf;  2.  immature  female;  3,  mature  female;  4.  mature  females  with 
newborn  calves  (mothers);  5,  dead;  6:  immature  male:  7,  mature  male.  Each  arrow 
represents  a  possible  transition  in  stage  from  one  year  to  the  next;  the  arrows  going  to 
stage  5  represent  stage- specific  mortalities.  A  calf  is  an  individual  that  was  sighted  along 
with  its  mother.  An  immature  is  an  individual  known  to  be  less  than  9  years  old.  Mature 
i  individuals  are  known  to  be  at  least  9  years  old,  or  in  the  case  of  females,  have  been 
j  sighted  previously  with  a  calf.  Mothers  are  females  that  are  sighted  with  a  newborn 
i  offspring.  If  the  sex  of  an  individual  is  unknown,  we  assume  it  has  an  equal  chance  of 
!  being  either  female  or  male,  as  our  data  contain  almost  equal  numbers  of  individuals  (1 41 
■  and  143)  known  to  be  female  and  male,  respectively.  Maturity  status  (whether  an 
individual  is  immature  or  mature)  is  unknown  in  22%  of  sightings.  When  maturity  status 
was  unknown,  we  estimated  the  probability  that  the  individual  was  immature  (0.30  for 
males,  0.87  for  temales)  using  the  method  described  in  ref.  10.  These  probabilities  were 
used  in  stage-assignment  matrices10  for  likelihood  calculations. 
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population8).  The  parameters  in  the  matrix  model  were  defined 
by  a  series  of  statistical  models,  the  parameters  of  which  were 
estimated  by  applying  stage-structured  mark-recapture  analysis5,910 
to  photographic  identification  data  collected  by  the  New  England 
Aquarium  (NEA)11.  The  best  model  was  selected  using  Akaike 
Information  Criteria  (AIC),2-M  (see  Methods). 


Figure  2  Stage-specific  sighting  and  survival  probabilities  for  males  and  females, 
a,  b.  Stage-specific  sighting  probabilities  of  the  best  model  (M?)  for  females  (a)  and  males 
(b).  c,  d.  Stage- specific  survival  probabilities  {M£  for  females  (e)  and  males  (d).  Squares, 
calf;  triangles,  immature;  circles,  mature;  diamonds,  mothers. 
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Figure  3  Constant  transition  probabilities  estimated  using  the  best  sighting  modal  («,)■ 
Life  cycle  graphs  ot  female  (a)  and  male  (b)  right  whales.  The  numbers  in  the  circles  are 
the  same  as  in  Fig.  1. 

In  the  resulting  model,  the  sighting  probability  of  most  stages 
depends  on  both  northern  and  southern  effort  levels  (Fig.  2a,  b;  see 
Methods  for  definition).  This  is  consistent  with  a  previous  study', 
which  found  a  significant  correlation  between  total  effort  and 
sighting  probability  in  a  simpler  model  that  did  not  distinguish 
sex  or  stages.  In  contrast  to  other  stages,  mothers  have  a  consistently 
high  sighting  probability  (Fig.  2a).  This  suggests  that  they  are  easier 
to  sight  than  other  stages,  and  that  the  survey  effort  in  the  two 
regions  has  been  sufficient  to  detect  almost  all  births  with  high 
probability. 

Using  the  best  sighting  probability  model,  we  estimated  two 
transition  models.  One  model  assumed  time-invariant  transition 
probabilities  (M,)  and  the  other  was  the  best  time-varying  model 
according  to  the  AIC  criterion  (M2).  The  time-invariant  model  (Mi, 
bp  =  0  for  all  j  and  i  in  equation  (2))  gives  a  weighted,  time-averaged 
picture  of  transition  probabilities.  In  this  model,  survival  prob¬ 
ability  is  low  for  calves,  higher  for  immature  and  mature  females 
and  much  lower  for  mothers  (Fig.  3).  According  to  the  best  time- 
varying  model  (M;),  the  transition  probabilities  of  mature  females 
and  males  have  been  constant,  whereas  the  other  stages  have 
transition  probabilities  that  change  with  time.  The  survival  prob¬ 
ability  of  mothers  shows  the  most  pronounced  decline  over  time 
(Fig.  2c,  d).  This  trend  is  statistically  significant,  as  the  slope 
parameter  (&„)  in  the  polychotomous  logistic  function  (equation 
(2))  is  significantly  below  0  (Fig.  4). 

A  series  of  population  projection  matrices,  A„  was  constructed  by 
augmenting  the  female  transition  probability  matrix  with  elements 
describing  reproduction"1  (see  Methods).  These  matrices  were  used 
to  calculate  population  growth  rate,  life  expectancy  and  expected 
lifetime  number  of  reproductive  events  experienced  by  a  female. 

The  asymptotic  population  growth  rate  (A)  is  given  by  the 
dominant  eigenvalue  of  the  population  projection  matrix.  When 
A  >  1,  the  population  is  asymptotically  growing;  when  A  <  1,  the 
population  is  asymptotically  declining.  Therefore,  A  is  an  important 
indicator  of  the  status  of  a  population.  If  we  assume  all  demographic 
parameters  have  been  constant  from  1980  to  1995  (M, ),  we  find  that 
A  =  1.01  (95%  confidence  interval  =  [1.00, 1.02);  the  confidence 
interval  was  estimated  using  the  estimated  covariance  of  parameters 
by  taking  the  inverse  of  the  Hessian  matrix"'),  suggesting  a  popula¬ 
tion  increase  of  1%  per  year,  on  average,  between  1980  and  1995. 
However,  time-specific  asymptotic  growth  rates  (A,) — calculated 
from  the  time- varying  matrices  A,  using  model  M2 — declined  from 
Alwo  =  1.03  to  Am,  =  0.98  (Fig.  5c).  Population  growth  rate 
declined  below  1  at  around  1992.  If  the  1995  population  growth 
rate  were  maintained,  the  population  would  go  extinct. 

The  decline  in  A,  is  the  net  result  of  all  the  changes  in  the  vital  rates 
(a  collective  term  for  transition,  survival,  fertility  and  mortality 
rates).  We  determined  how  much  the  decline  in  each  vital  rate 
contributed  to  the  decline  in  A„  using  a  life  table  response  experi¬ 
ment  (LTRE)  analysis5*.  This  analysis  decomposes  the  changes  in  A, 
into  contributions  from  each  entry  in  A,.  Figure  5d  shows  the  sum  of 
the  contributions  of  all  matrix  entries  involving  each  stage.  The 


Figure  4  Joint  95*  profile  likelihood  confidence  region  lor  the  logistic  parameters  of 
survival  probability  ot  mothers. 

results  of  this  analysis  show  that  the  decline  in  A,  is  caused  mainly  by 
the  decline  in  the  survival  probability  of  mothers  (Fig.  5d). 

The  life  expectancy  at  birth  is  the  mean  age  at  which  individuals 
bom  in  a  given  year  would  die  if  conditions  of  that  year  were 
maintained.  It  can  be  calculated  by  treating  the  transition  matrix  as 
an  absorbing  Markov  chain  and  calculating  the  mean  time  to 
absorption  (that  is,  death)5.  During  the  early  1980s,  the  life 
expectancy  of  females  was  twice  that  of  males  (Fig.  5a).  This  may 
be  typical  of  large  cetaceans.  (Killer  whales  (Orcinus  orca),  for 
example,  exhibit  a  similar  difference  in  life  expectancy  between 
females  and  males15.)  Female  life  expectancy  has  declined  from 
about  51.8  years  in  1980  to  about  14.5  years  in  1995.  Until  recently, 
life  expectancy  of  females  has  exceeded  that  of  males,  but  that  is  no 
longer  true  owing  to  the  reduced  survival  probability  of  females. 

The  expected  number  of  reproductive  events  during  a  female’s 
lifetime5  has  declined  from  about  5.27  in  1980  to  about  1.26  in  1995 
(Fig.  5b).  A  mature  female  could  once  expect  to  reproduce  6.48 
times;  that  number  is  now  about  1.80. 

The  growth  rate  projections  in  Fig.  5c  are  deterministic.  The  right 
whale  population  is  small,  however,  and  the  population  of  any  single 
stage  is  even  smaller.  Because  of  this,  it  has  been  suggested  that 
population  projections  should  include  demographic  stochasticity 
(chance  fluctuations  due  to  the  random  fates  of  a  small  number  of 
individuals;  see  ref.  5).  To  evaluate  the  effect  of  demographic 
stochasticity,  we  transformed  the  matrix  A|TO  into  a  multitype 
branching  process5,  and  calculated  the  probability  distribution  of 
times  to  extinction  (Fig.  6).  We  compared  this  distribution  with  the 
deterministic  prediction  obtained  by  using  and  defining 

extinction  as  the  lime  when  total  population  size  reached  1 
individual.  For  both  calculations,  we  used  an  initial  condition  of 
150  females  distributed  according  to  the  stable  stage  distribution. 
This  simulates  the  hypothetical  situation  in  which  the  vital  rates  of 
1995  remain  constant  and  neither  environmental  trends  nor  envir¬ 
onmental  fluctuations  affect  them. 

The  mean  time  to  extinction  under  demographic  stochasticity  is 
208  years  (similar  to  the  estimate  of  191  years  in  ref.  1).  The 
deterministic  lime  to  extinction  is  245  years.  Thus  demographic 
stochasticity  reduces  expected  extinction  time  by  15%.  Figure  6 
gives  the  complete  distribution  of  extinction  times;  there  is  a  5% 
chance  of  extinction  within  130  years  and  a  25%  chance  of  extinc¬ 
tion  within  165  years.  These  calculations  exclude  other  factors  such 
as  continued  declining  survival  trends,  environmental  stochasticity 
and  Allee  effects,  all  of  which  would  hasten  extinction.  Thus  the 
expected  time  of  208  years  should  be  considered  an  upper  bound. 

Right  whale  conservation  efforts  are  directed  towards  reducing 
mortality  due  to  entanglement  and  ship  collisions.  Because  the 
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Figure  5  Demographic  parameters  calculated  from  the  time-varying  matrices  A,,  a.  The 
mean  life  expectancy  at  birth  for  males  (triangles)  and  females  (circles),  b.  The  mean 
number  of  lifetime  reproductive  episodes,  estimated  from  birth  (circles)  and  from  maturity 
(triangles),  c,  Asymptotic  population  growth  rates  \f  (solid  line)  calculated  from  the 
population  projection  matrices  A,  produced  by  the  best  model  (Ay  and  from  the 
unstructured  model  of  ref.  1  (dotted  line).  The  error  bars  were  standard  errors,  which  were 
calculated  using  the  series  approximation  to  the  variance  of  X  and  the  covariance  matrix 
obtained  from  the  information  matrix,  d.  The  result  of  a  LTRE  decomposition  analysis  for 
the  trend  in  X.  The  demographic  rates  in  1980  are  used  as  the  reference  matrix,  and  the 
contributions  of  all  matrix  entries  involving  each  stage  were  summed.  The  solid  line  shows 
the  actual  trend  in  X,  measured  as  a  deviation  from  the  value  in  1980.  It  is  closely 
approximated  by  the  sum  of  all  the  other  lines.  Squares,  calf;  triangles,  immature;  circles, 
mature;  diamonds,  mothers. 
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Extinction  time  (year) 

Figure  6  The  probability  distribution  of  time  to  extinction  assuming  demographic 
stochasticity.  The  distribution  was  calculated  from  a  multi-type  branching  process, 
treating  transitions  and  reproduction  as  independent  events  (multinomial  and  Bernoulli, 
respectively;  see  ret.  5). 

population  is  so  small,  a  single  death  represents  a  significant 
mortality  rate.  Conversely,  significant  reduction  in  mortality  rate  = 
can  be  obtained  by  preventing  just  a  few  deaths.  Figure  7  shows  the 
effect  on  X  of  preventing  female  deaths,  using  the  vital  rates  of  1995 
and  assuming  a  starting  population  of  150  females  distributed 
according  to  the  stable  stage  distribution.  The  results  indicate  that 
prevention  of  just  two  female  deaths  every  year  would  suffice  to 
increase  X  to  1.  This  gives  a  sense  of  the  magnitude  of  the  initial 
management  action  needed  for  the  protection  of  the  population. 
However,  as  population  size  increases  owing  to  successful  manage¬ 
ment,  more  deaths  will  have  to  be  prevented  to  maintain  a  positive 
population  growth  rate  because  the  number  of  deaths  that  translates 
into  a  given  mortality  rate  also  increases  with  the  population  size. 

The  causes  behind  the  decline  in  survival  probability  of  mothers 
are  still  unknown;  however,  collisions  with  ships,  entanglement  with 
fishing  gear  and  changes  in  food  availability  due  to  climate  fluctua¬ 
tions  are  suspected  to  contribute  towards  mortality16,17.  Although 
right  whales  are  distributed  from  the  coast  of  northern  Florida  to 
the  Bay  of  Fundy,  it  is  primarily  females  and  calves  that  are  sighted 
in  the  calving  ground  off  the  coast  of  Florida  and  Georgia1".  The 
calving  ground  is  close  to  shipping  lanes  where  large  vessel  traffic 
has  increased  since  1980  (ref.  19).  More  than  60%  of  North  Atlantic 
right  whales  have  scars  from  entanglement  in  fishing  gear  such  as 
lobster  pots  and  sink  gillnets17.  We  have  found  a  significant  rank 
correlation  between  crude  survival  probability1  (the  survival  prob¬ 
ability  of  individuals  without  distinguishing  stage  differences)  and 
North  Atlantic  Oscillation  index  (NAO),  which  is  often  used  as  an 
index  for  climate  fluctuations20.  Studies  that  focus  on  the  effects  of 
these  three  factors  are  urgently  needed. 

This  study  combines  multi-stage  mark- recapture  methods  and 
matrix  population  model  analyses.  Both  methods  incorporate 
differences  in  vital  rates  that  are  experienced  by  different  stages 
within  a  population.  Multi-stage  mark-recapture  analysis  also 
incorporates  stage-specific  heterogeneity  in  sighting  probability  : 
within  the  population.  We  succeeded  in  identifying  a  life  cycle  ; 
stage  that  is  experiencing  reduced  survival  probability,  and  were 
able  to  use  the  model  to  document  the  implications  of  that  reduced 
survival  probability  for  population  growth. 

Like  any  mark- recapture  analysis,  ours  can  be  influenced  by 
heterogeneity  in  sighting  probability  among  individuals  within  a 
stage.  If  the  sighting  probability  of  mature  females  (stage  3)  was 
heterogeneous,  the  survival  probability  of  mothers  towards  the  end 
of  the  study  period  would  be  underestimated.  However,  hetero¬ 
geneity  in  sampling  alone  cannot  explain  the  observed  decline  in  the 
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Figure  7  The  predicted  population  growth  rate  that  would  result  from  preventing  deaths  of 
females  regardless  of  their  stage. 


survival  probability.  Furthermore,  during  the  winter  of  1996,  an 
unusually  large  number  of  confirmed  mortalities  were  reported  in 
the  southeast  United  States,  which  is  the  only  known  calving 
ground17.  This  evidence  suggests  that  the  declining  trend  in  the 
survival  of  mothers  is  both  real  and  a  great  concern. 

The  North  Atlantic  right  whale  is  currently  seriously  endangered 
owing  to  declining  survival  probability,  especially  among  mothers. 
This  finding,  together  with  the  unprecedented  calf  production  in 
the  spring  of  2001  (ref.  21).  suggests  that  Allee  effects  are  not  yet 
limiting  this  population.  Our  analysis  also  shows  that  the  popula¬ 
tion  was  experiencing  a  positive  growth  rate  in  the  early  1980s 
(Fig.  5c).  This  implies  that  it  is  not  necessary  to  return  the  vital  rates 
to  those  of  the  pre-whaling  period  to  obtain  a  positive  growth  rate. 
There  is  every  reason  to  hope  that  prompt  management  interven¬ 
tion  can  improve  survival  enough  to  permit  the  recovery  of  the 
North  Atlantic  right  whale.  □ 

Methods 

Photographic  identification  data 

Since  19B0,  iIk  NEA  ha*  accumulated  photographs  of  approximately  350  right  whales 
taken  on  more  than  10,000  sighting  occasions.  Individual  right  whales  can  be  recognized 
by  natural  markings  such  as  scars  and  callosity  patterns.  On  the  basis  of  these 
identifications  the  NEA  has  constructed  a  database  of  individual  sighting  histories,  which 
consist  of  information  on  whether  or  not  each  individual  was  sighted  at  least  once  in  each 
year.  These  data  can  be  treated  as  if  individuals  were  marked  on  the  occasion  of  their  first 
identifications,  and  recaptured  at  subsequent  sightings.  We  used  the  data  from  1980  to 
1996,  which  were  complete  at  the  beginning  of  this  analysis. 


,  Statistical  model* 

Sighting  probabilities  of  each  stage  and  transition  probabilities  among  stages  were 
estimated  by  fitting  a  series  of  statistical  models  to  the  sighting  history  data  by  maximum 
likelihood.  The  statistical  models  describe  sighting  probability  as  a  function  of  sampling 
effort  in  two  regions.  The  northern  effort  level  is  measured  by  the  total  annual  survey  days 
in  Massachusetts  Bay,  Great  South  Channel,  Bay  of  Fundy,  and  in  Brown’s  Bank;  the 
southern  effort  level  is  measured  by  the  total  annual  survey  days  off  the  coast  of  Florida 
and  Georgia.  The  dependence  of  sighting  probability  on  the  northern  effort,  the  southern 
effort,  both  efforts,  or  neither  effort  was  modelled  with  binary  logistic  functions.  Let  Si(r) 
;  be  the  sighting  probability  of  stage  i  at  time  t.  Then 


cxp(x.  +  y,tu  -t-  x,e:j) 

I  +  expfx,  +  y.e,.  +  z.e,,) 


(J> 


where  x,  is  an  intercept  parameter,  and  y,  and  z,  are  slope  parameters  associated  with 
northern  (e,.,)  and  southern  (e^)  effort  levels,  respectively. 

We  modelled  the  transition  probabilities  of  each  stage  as  poJychotomous  logistic 
functions  of  time”0'.  Let  p„(f)  be  the  transition  probability  from  stage  /  at  time  t  to; 
at  time  t+  I.  Then 


exp (a,  +  b,t) 


(2) 


where  a„  and  bf  are  intercept  and  slope  parameters,  respectively.  When  all  the  slope 
parameters  are  set  to  zero,  the  transition  probabilities  arc  time-invariant. 

Model  selection 

Model  selection  was  done  in  the  following  sequence.  First.  1.024  sighting  models  were 
created  by  including  all  possible  combinations  of  effort  levels  for  all  possible  combinations 
of  stages  (equation  (1)).  These  models  were  combined  with  the  transition  model  in 
equation  (2)  m  which  all  probabilities  arc  functions  of  lime.  We  selected  the  best  of  these 
models  using  A1C,:.  Tbc  A1C  difference  between  the  best  and  the  second-best  sighting 
models  was  about  2,  indicating  that  the  support  for  the  best  model  relative  to  the  second 
best  model  is  high1*.  Therefore,  we  used  only  the  best  model.  Because  sighting  probabilities 
of  immature  males  and  females  did  not  differ  significantly,  by  a  likelihood  ratio  test,  they 
were  set  equal.  Using  the  resulting  model  for  sighting,  we  selected  the  best  time-varying 
transition  matrix  from  all  64  models  created  by  allowing  the  transition  probabilities  of  all 
possible  combination  of  stages  to  depend  on  time  according  to  equation  (2).  The  four  best 
transition  models  had  AIC  differences  from  the  best  model  of  less  than  2,  indicating  that 
the  data  provide  some  support  for  all  these  models.  All  of  these  models,  however,  had 
time-dependent  survival  probability  of  mothers. 

Projection  mafrlx 

The  projection  matrix  is 


where 


(0  r,  0  \ 

p„l’>  0  0  I 

0  Pull)  P»l‘)  Pi.H>  I 

(I  p,._  It)  PulD  0  I 

r.  =  o Sp„u)p’J(t  +  D 


and 


F,  =  0.5/>41<r)p?;<r  +  1) 


(3) 


(4) 


(S) 


The  first  row  of  A,  describes  reproduction;  the  other  entries  are  transition  probabilities. 
Consider  F2.  When  a  female  moves  from  stage  2  to  stage  4  (with  probability  />4,<M).  she 
gives  birth;  the  newborn  calf  is  female  with  probability  0.5.  Although  newborn  calves  have 
distinct  markings,  they  are  harder  to  distinguish  individually  than  other  stages.  There¬ 
fore.  calf  survival  is  estimated  from  the  point  where  the  calf  is  seen  sufficiently  well  to 
permit  identification,  which  is  not  necessarily  on  its  first  sighting.  To  appear  as  a  calf  in 
stage  1  at  I  +  I,  the  newborn  calf  must  survive  long  enough  to  be  catalogued.  We  assume 
calves  are  catalogued  on  average  midway  through  their  first  year,  and  that  the  mother 
must  survive  this  long  (with  probability  pVl r  +  I ))  in  order  for  the  calf  to  survive.  F»  is 
similar. 
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Chapter  6 

Stage-structured  demographic 
analysis  of  North  Atlantic  right 
whale:  Effects  of  environmental 
and  demographic  stochasticity  on 
population  viability 

6.1  Introduction 

North  Atlantic  right  whales  ( Eubalaena  glacialis )  are  one  of  the  most  endangered  mammal 
populations  in  the  world.  By  the  beginning  of  the  20 th  century,  the  population  in  the 
western  Atlantic  had  been  hunted  to  near  extinction.  Despite  more  than  a  half  century  of 
international  protection,  no  clear  sign  of  their  recovery  has  been  detected.  This  population 
currently  consists  of  fewer  than  300  individuals,  and  there  is  evidence  to  suggest  that  this 
number  is  slowly  shrinking  because  of  declining  survival  probability  (Caswell  et  al.,  1999, 
also  Chapter  5).  A  previous  study  (Chapter  5)  has  shown  that  this  reduction  in  survival  is 
primarily  a  result  of  an  increase  in  the  mortality  of  females  during  the  year  after  giving  birth. 
The  same  study  has  also  concluded  that  if  the  environmental  conditions  experienced  by  this 
population  in  1995  were  maintained  in  future  years,  the  population  would  be  expected  to  go 
extinct  in  about  200  years.  The  survival  of  this  population  depends  on  prompt  and  effective 
conservation  efforts  based  on  our  understanding  of  the  population. 
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Results  of  any  demographic  analysis  has  inherent  uncertainties  because  a  population 
is  affected  by  many  factors;  the  North  Atlantic  right  whales  are  not  an  exception.  For 
example,  demographic  stochasticity  can  reduce  the  time  to  extinction  (Caswell,  2001).  The 
results  in  Chapter  5  and  Pershing  &  Greene  (2002)  showed  that  the  estimated  survival 
probability  and  reproductive  rate,  respectively,  are  significantly  correlated  with  an  index 
for  climate  pattern  in  the  North  Atlantic  (the  North  Atlantic  Oscillation  index),  suggesting 
that  the  environmental  forcing  is  also  affecting  the  population.  In  this  paper,  we  present 
the  results  of  a  demographic  analysis  of  the  North  Atlantic  right  whale  that  analyzes  effects 
of  demographic  stochasticity  and  environmental  stochasticity  on  population  viability. 

Another  uncertainty  inherent  in  demographic  analyses  are  the  uncertainties  in  parameter 
estimation  and  model  selection.  In  this  paper,  we  will  assess  the  parameter  estimation 
uncertainty  using  a  parametric  bootstrap  method.  We  also  apply  the  Akaike  Information 
Criteria  (AIC)  to  select  the  best-fit  models  from  multiple  candidate  models.  This  method 
takes  into  account  model  selection  uncertainty. 

6.2  Method 

6.2.1  Matrix  Population  Model  for  the  North  Atlantic  right  whales 

Our  analysis  uses  a  stage-structured  matrix  population  model  (e.g.  Caswell,  2001).  This 
model  takes  into  account  the  differences  in  vital  rates  among  life  stages.  The  model  in  a 
general  form  is  expressed  as 

n(t  +  1)  =  A  ^n(f)  (6.1) 

where  n(t )  is  a  vector  containing  the  number  of  individuals  at  year  t,  and  A-^  is  a  matrix 
that  projects  the  number  of  individuals  in  each  stage  from  year  t  to  t  +  1.  The  entries  in 
A^  are  determined  by  a  life  cycle  structure  of  a  population  of  interest. 

For  the  North  Atlantic  right  whales,  we  developed  a  life  cycle  structure  based  on  available 
information  on  their  life  history  Figure  6.1.  In  the  figure,  the  numbers  indicate  stages,  solid 
arrows  indicate  possible  transitions  that  an  individual  can  make  from  one  year  to  the  next 
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Figure  6.1:  Life  cycle  graph  of  the  North  Atlantic  right  whales:  1:  calves,  2:  immature 
females,  3:  mature  females,  4:  mothers,  5:  mature  females  during  a  year  after  the  mother 
stage 


if  it  survives,  and  the  dashed  arrow  represents  reproduction  of  offspring  by  mature  females 


(i.e.  fertility).  Each  of  the  solid  arrows  has  an  associated  transition  probability  4j  ,  which 
gives  the  probability  that  an  individual  in  stage  j  in  year  t  is  in  stage  i  in  the  following 
year.  This  stage  structure  is  slightly  different  from  the  one  in  Chapter  5.  In  the  new  stage 
structure,  females  have  an  extra  stage  after  giving  birth  (stage  5)  during  which  they  cannot 
reproduce  again.  This  reflects  the  fact  that  individual  right  whales  rarely  give  birth  within 
two  years  after  a  successful  reproduction. 

Based  on  the  stage  structure  in  Figure  6.1,  we  wrote  down  a  projection  matrix  for  the 
North  Atlantic  right  whales: 
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(6.2) 


The  first  row  of  describes  reproduction;  the  other  entries  are  the  transition  probabilities. 
Here,  we  assume  that  population  dynamics  are  determined  only  by  female  vital  rates  and 
there  are  always  enough  males  to  fertilize  them.  Therefore,  we  only  keep  track  of  number 
of  females. 
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Because  of  the  way  the  stage  structure  is  defined,  the  fertility  Fl  ^  can  be  expressed  in 
terms  of  transition  probabilities  as 

=  0.5 (4+1))°'5  •  (6-3) 

This  relationship  is  based  on  the  following  arguments.  When  a  female  moves  from  stage  3 
to  stage  4  (with  probability  <j>\ 2),  she  gives  birth;  the  calf  is  female  with  probability  0.5. 
Although  calves  have  distinct  markings,  they  are  more  difficult  to  distinguish  individually 
than  other  stages.  Therefore,  calf  survival  is  estimated  from  the  point  where  the  calf  is  seen 
sufficiently  well  to  permit  identification,  which  is  not  necessarily  on  its  first  sighting.  To 
appear  as  a  calf  in  stage  1  in  year  f  +  1,  the  newborn  calf  must  survive  long  enough  to  be 
cataloged.  We  assume  calves  are  cataloged  on  average  midway  through  their  first  year,  and 
that  the  survival  probability  is  equal  to  that  of  mother  during  the  same  period  (with  prob¬ 
ability  (<jf>54+1))°  5)-  Ordinarily,  when  estimating  a  population  projection  matrix,  transition 
probabilities  and  fertility  terms  must  be  estimated  separately,  often  requiring  different  sets 
of  data.  Because  we  can  express  the  fertility  term  with  the  transition  probabilities,  we  only 
need  data  to  estimate  the  transition  probabilities. 

6.2.2  Data 

The  data  we  used  to  estimate  transition  probabilities  are  individual  sighting  histories  of 
the  right  whales.  They  are  based  on  photographs  taken  by  members  of  the  North  Atlantic 
right  whale  consortium  along  the  east  coast  of  the  United  States  from  1980  to  1998.  Major 
survey  areas  include  Massachusetts  Bay,  Great  South  Channel,  Bay  of  Fundy,  Brown’s  Bank, 
and  off  the  coast  of  Florida  and  Georgia.  Because  right  whales  have  unique  markings  called 
callosity  patterns  on  their  head  and  scars  from  entanglement  in  fishing  gear,  individuals  can 
be  identified  from  the  photographs.  Based  on  these  photographs,  New  England  Aquarium 
has  been  accumulating  annual  sighting  histories  of  individual  right  whales.  For  the  purpose 
of  our  analysis,  individuals  were  considered  captured  and  marked  at  the  occasion  of  their 
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first  sighting  and  recaptured  when  they  were  resighted  at  least  once  in  a  later  year. 

Of  the  406  individuals  captured  between  1980  and  1998,  158  were  identified  as  female 
and  163  were  identified  as  male.  The  rest  were  assumed  to  be  female  or  male  with  an  equal 
probability.  At  each  sighting  occasion,  we  also  categorized  the  sighted  individual  into  one 
of  the  five  life  stages  (calf,  immature,  mature,  mother,  and  mature  female  during  a  year 
after  the  mother  stage)  based  on  the  following  criteria.  Calves  were  sighted  along  with  their 
mothers.  Immature  individuals  were  known  to  be  less  than  9  years  old.  Mature  individuals 
were  known  to  be  at  least  9  years  old,  or  in  the  case  of  females,  had  been  sighted  previously 
with  a  calf.  A  “mother”  is  a  female  during  the  year  after  a  successful  birth.  A  stage  for 
mature  female  during  a  year  after  the  mother  is  for  a  female  that  is  known  to  have  given 
a  birth  in  a  previous  year.  After  assigning  the  stages  based  on  these  criteria,  we  were  left 
with  about  30%  of  the  total  sightings  that  belong  to  either  immature  or  mature  stages  but 
did  not  have  enough  information  to  assign  either  stage  with  certainty.  We  assumed  that 
they  have  an  equal  chance  of  being  immature  or  mature. 

6.2.3  Multi-stage  mark-recapture  statistics 

We  estimate  the  transition  probabilities  <f>fj  from  the  sighting  (capture)  history  sequences. 
A  difficulty  in  estimating  transition  probabilities  from  the  capture  history  sequences  such 
as  the  ones  available  for  the  North  Atlantic  right  whales  comes  from  the  fact  that  not  all 
individuals,  even  if  they  are  still  alive,  are  captured  every  year.  If  a  capture  history  sequence 
ends  without  capturing  an  individual  for  one  or  more  sampling  occasions,  we  do  not  know 
whether  the  individual  was  dead  or  alive.  It  could  have  died  immediately  after  the  last 
capture,  it  could  have  died  later,  or  it  may  still  be  alive.  Furthermore,  if  an  individual 
was  not  captured,  we  often  do  not  have  enough  information  to  determine  the  stage  of  that 
individual  for  sure. 

Instead  of  assuming  any  particular  whale  is  dead  or  alive,  we  modeled  the  capture  history 
sequence  using  both  capture  probabilities  pf  \  which  give  the  probability  of  capturing  an 
individual  that  is  in  stage  i  in  year  t,  and  transition  probabilities  4>if  ■  With  the  capture  and 
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transition  probabilities,  we  calculated  the  sum  of  the  probabilities  of  every  possible  sequence 
of  events  that  was  compatible  with  a  given  capture  history  to  obtain  the  probability  of  that 
capture  history.  The  product  of  all  the  capture  history  probabilities  was  reinterpreted  as 
a  likelihood.  Finally  a  maximum  likelihood  method  was  applied  to  estimate  the  transition 
probabilities  as  well  as  the  capture  probabilities  (see  Nichols  et  al.,  1992b,  also  Chapter  2). 

Major  assumptions  of  the  multi-stage  mark-recapture  method  in  general  include  (1) 
individuals  within  the  same  stage  make  transitions  into  another  stage  independently  but 
based  on  the  same  probability  and  (2)  marked  individuals  within  the  same  stage  are  captured 
independently  but  based  on  the  same  probability.  However,  we  do  not  assume  that  any 
particular  whale  is  dead  or  alive,  nor  do  we  assume  that  unmarked  individuals  have  the 
same  capture  probability  as  marked  individuals.  For  the  right  whale  study  in  particular,  we 
also  assume  that  (1)  when  the  stage  of  an  individual  is  unknown,  it  is  immature  or  mature 
with  equal  probability,  (2)  when  the  sex  of  an  individual  is  unknown,  it  is  female  or  male 
with  equal  probability,  and  (3)  the  capture  probability  of  a  dead  individual  is  0.  These 
assumptions  are  incorporated  into  the  likelihood  using  the  method  described  in  Chapter  2. 

6.2.4  Parameters  to  be  estimated 

Figure  6.2  shows  a  complete  transition  graph  of  the  North  Atlantic  right  whale.  This 
transition  graph  includes  male  stages  and  another  stage  for  death.  Stages  in  this  transition 
graph  correspond  to  the  stages  we  assigned  to  the  sighted  individuals.  Each  arrow  indicates 
all  possible  stage  transitions  (including  death)  of  individuals  from  one  year  to  the  next  and 
has  associated  transition  probability  Once  we  estimate  the  transition  and  capture 

probabilities  associated  with  this  transition  graph,  we  can  use  the  transition  probabilities 
associated  with  female  stages  to  construct  the  population  projection  matrix  (Equation  6.2). 
We  let  the  survival  probability  of  calves  between  females  and  males  to  be  the  same,  and  we 
need  male  data  to  estimate  this  probability.  Because  of  this  assumption,  we  still  need  to 
estimate  parameters  associated  with  males. 

From  the  19  year  capture  history  sequences,  we  estimated  parameters  that  define  transi- 
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tion  and  capture  probabilities.  The  transition  graph  in  Figure  6.2  contain  8  live  stages  and 
19  possible  stage  transitions.  If  all  probabilities  were  free  to  vary  from  one  year  to  the  next, 
there  would  be  about  300  parameters  to  estimate.  The  number  of  parameters,  however, 
can  be  reduced  by  making  reasonable  assumptions  about  how  the  probabilities  are  related 
among  stages  and  on  how  they  change  with  time. 


Figure  6.2:  Stage  transition  graph  of  the  North  Atlantic  right  whales:  0:  deaths,  1:  calves, 
2:  immature  females,  3:  mature  females,  4:  mothers,  5:  mature  females  during  a  year  after 
the  mother  stage,  6:  immature  males,  7:  mature  males 
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Capture  probability  model 


To  reduce  the  number  of  parameters  associated  with  capture  probabilities,  we  assumed 
that  they  change  freely  from  one  year  to  the  next  but  are  correlated  between  stages.  Let 
(pW  be  the  logistic  transform  of  the  capture  probability  of  mature  females  (stage  3)  (i.e. 
<pW  =  log  We  let  logistic  transforms  of  capture  probabilities  of  other  stages  to 

be  linear  functions  of  tp®  (i.e.  ipi  +  Ui<pW  =  log  ~^y).  After  we  transform  (p®  and  the 
functions  back  to  probabilities,  we  obtain  the  following  statistical  model: 


Ps(<) 

Pi(t ) 


exp  (<pW) 

1  -f  exp  (<pW)  ’ 
exp  (ipj  + 

1  +  exp  (ipi  +  l 


for  i  ^  1  or  3, 


(6.4) 

(6.5) 


where  t  is  year  (for  capture  probability  t  =  1981, . . . ,  1998),  and  </?*,  ipi,  and  u >i  are  parame¬ 
ters  to  be  estimated.  Equation  (6.4)  is  a  logit  function  with  a  different  parameter  for  each 
year.  Equation  (6.5)  is  a  logistic  function  with  ip*  as  a  covariate.  In  this  equation,  ipi  is  an 
intercept  and  u>i  is  a  slope  parameter.  Because  the  capture  probabilities  other  than  mature 
females  are  expressed  only  with  two  parameters  ipi  and  cuj,  this  model  reduces  the  number 
of  parameters. 

We  used  stage  3  (mature  female)  as  a  reference  stage,  for  which  one  parameter  is  es¬ 
timated  for  each  year  (see  equation  (6.4)),  because  the  data  contain  this  stage  the  most. 
The  rest  of  stages  use  the  parameters  in  equation  (6.4)  as  a  covariate.  This  model  does  not 
include  the  capture  probability  of  calves.  It  is  not  needed  for  estimating  transition  proba¬ 
bilities  because  mark-recapture  statistics  are  conditional  on  marking  of  animals.  Capturing 
of  calves  always  occurs  at  the  very  first  capture  (i.e.  marking)  of  the  individuals. 

The  resulting  capture  probability  model  contains  28  parameters  to  be  estimated.  We 
used  this  model  for  all  of  the  analyses  presented  in  this  paper. 
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Transition  probability  models 


In  the  current  study,  we  focus  on  modeling  the  reproductive  rate  and  the  stage  specific 
survival  probability  of  calves  and  mothers.  A  previous  study  in  Chapter  5  has  shown  that 
the  survival  probability  of  mothers  and  calves  has  declined  significantly  over  time;  on  the 
other  hand,  the  survival  probability  of  other  stages  remained  constant.  The  same  study 
has  also  concluded  that  the  survival  probability  of  individuals  averaged  over  all  stages  in 
Caswell  et  al.  (1999)  is  significantly  correlated  with  the  NAO  index  that  is  lagged  by  2 
years.  Furthermore,  Pershing  &  Greene  (2002)  have  shown  that  reproductive  rate  of  the 
population  has  been  significantly  correlated  with  the  same  NAO  index.  Therefore,  we  model 
transition  probabilities  incorporating  the  trend  associated  with  time  and  the  NAO  index. 

In  order  to  model  the  reproductive  rate  and  the  survival  probability  of  mature  fe¬ 
males  separately,  we  decomposed  the  transition  probability  </> 43  into  a  stage-specific  sur¬ 
vival  probability  and  a  transition  probability  conditional  on  the  survival  of  individuals  (i.e. 
<42  =  4^42)-  We  call  the  latter  a  conditional  transition  probability,  and  these  two  proba¬ 
bilities  are  denoted  by  s ^  and  42  >  respectively.  The  stage-specific  survival  probability  gives 
the  probability  of  individuals  in  stage  3  at  year  t  surviving  to  the  following  year,  and  the 
conditional  transition  probability  gives  the  probability  that  individuals  in  stage  3  in  year  t 
will  be  in  stage  4  in  the  following  year  conditional  on  the  individuals  survival  through  the 
year.  The  relationship  between  A3  and  these  two  probabilities  are 

AS  =  i-44).  (6-g) 

<42  -  40u-42)>  M 

<42  -  4X1  (e-s) 

In  the  above  formulation,  42  *s  a  probability  of  mature  females  giving  birth  conditional  on 
their  survival.  Thus  it  gives  the  reproductive  rate  of  mature  females. 

The  reproductive  rate  (42)  and  the  survival  probabilities  of  calves  (s^  =  <4i)  and 


mothers  (s^  =  <$})  are  modeled  with  logistic  functions  of  two  covariates: 


»i(*)  = 
s4(t)  = 

f43(i)  = 


exp  (qi  -f  Pi cf1  +  7ic^f)) 

1  +  exp  (ai  +  A  4^  +  7W2  4 
exp  (a2  +  A4f)  +  I1C2) 

1  +  exp  (a2  +  +  72C2  4 

exp  (a3  +  Ascj^  +  73C2 ' >) 

1  +  exp  (a3  +  +  73  c2  4 


(6.9) 

(6.10) 

(6.11) 


where  t  is  year  (for  transition  probabilities  t  =  1980, . . . ,  1997),  and  are  values  of 
the  winter  NAO  index  (December  of  year  t-  3  to  March  of  year  t-  2)  and  a  time  index 
(0,  Yf,  ^ , . . . ,  , . . . ,  1)  for  year  t,  ak  (k  =  1, 2, 3)  is  an  intercept  parameter,  and  (3k  and 

7^  are  slope  parameters  associated  with  the  covariates.  By  setting  some  or  all  of  (3k  and  7 k 
to  zeros,  we  constructed  64  models  that  differ  in  the  dependency  of  the  probabilities  on  the 
two  covariates. 

To  be  consistent  with  the  transition  probability  of  mature  females,  we  also  decompose 
the  transition  probabilities  of  immature  stages  (stage  2  and  6)  into  stage-specific  survival 
probabilities  (s2  and  sg)  and  conditional  transition  probabilities  (t32  and  t^).  However, 
these  probabilities  along  with  the  survival  probabilities  of  mature  females  and  males  are 
assumed  to  be  constant  over  time.  For  each  constant  probability,  its  logistic  transform  was 
estimated  as  a  parameter. 

We  fitted  the  64  models  to  the  capture  history  sequences,  calculated  AIC  for  each  model, 
and  compared  the  AIC  values  to  select  the  best  fit  model. 


Saturated  model 

In  addition  to  the  64  models,  we  also  used  another  model  in  which  all  transition  probabilities 
vary  freely  from  one  year  to  the  next.  We  call  this  model  a  saturated  model.  This  model 
contains  a  total  of  184  parameters  (156  are  associated  with  transition  probabilities,  27  are 
associated  with  capture  probabilities,  and  1  is  a  confounded  parameter  that  is  a  product 
of  stage-specific  survival  probability  from  1997  to  1998  and  capture  probability  of  mature 


105 


females  in  1998). 

We  assume  that  any  deviations  from  the  best  models  among  the  64  models  are  resulting 
from  effects  of  environmental  variability.  And  we  used  the  saturated  model  to  calculate 
deviations  in  parameter  values  due  to  the  environmental  stochasticity. 


6.3  Results 

6.3.1  Model  selection 

Prom  the  64  models,  we  select  the  best  model  that  optimize  the  errors  of  over-fitting  by 
including  too  many  parameters  and  under-fitting  by  over-simplifying  models.  This  model 
selection  was  done  using  Akaike  information  criteria  (AIC)  (Burnham  &  Anderson,  1998). 
AIC  is  an  index  for  weight  of  evidence  in  the  data  that  a  corresponding  model  is  the  best 
model  among  the  ones  compared.  AIC  values  associated  with  64  models  were  transformed 
into  Akaike  weights,  which  are  the  normalized  AIC  indices  aso  that  they  sum  to  1.  Using 
the  index,  we  chose  the  10  best  models  (Table  6.3.1).  In  the  same  table,  I  also  show  the 
difference  between  AIC  for  each  model  and  one  associated  with  the  best  model  (i.e.  A 
AIC).  The  best  and  the  second  best  models  are  supported  almost  equally  at  0.27  and  0.25, 


Table  6.1:  Akaike  weight  for  the  10  best  models,  y/  indicates  that  the  model  include  the 
corresponding  covariate  (time  or  NAO  index). 


Model 

Akaike 

A  AIC 

Calf  Survival 

Mother  Survival 

Reproduction 

Weight 

NAO 

Time 

NAO 

Time 

NAO 

Time 

WSM 

0 

V 

7 

7 

7 

US 

0.18 

V 

V 

7 

7 

m3 

US 

1.24 

V 

7 

y/ 

7 

7 

Mi 

0.10 

2.06 

V 

7 

y/ 

s/ 

7 

m5 

0.07 

2.78 

7 

7 

7 

7 

7 

M6 

0.05 

3.47 

V 

7 

V 

yj 

7 

7 

M7 

0.03 

4.50 

V 

7 

7 

M8 

0.02 

5.29 

7 

7 

m9 

0.02 

5.30 

7 

7 

7 

Mio 

0.01 

6.40 

7 

7 

7 

7 
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respectively.  In  both  models,  the  survival  probability  of  mothers  is  correlated  with  NAO 
and  the  reproductive  rate  is  correlated  with  both  the  time  and  NAO  indices.  However, 
these  models  differ  in  how  survival  probability  of  calves  changes  over  time.  The  best  model 
suggests  that  the  survival  probability  depends  on  the  time  index  and  is  monotonically 
declining,  but  the  second  model  predicts  that  the  survival  probability  is  correlated  with  the 
NAO  index. 

The  AAIC  values  for  the  first  six  models  in  Table  6.3.1  are  less  than  4  suggesting 
these  models  deserve  our  attention  (Burnham  &  Anderson,  1998).  In  all  six  models,  the 
reproductive  rate  is  correlated  with  the  NAO  index  indicating  that  there  is  a  strong  evidence 
in  the  data  that  the  environmental  forcing  is  influencing  right  whale  reproduction.  In  the  six 
models,  the  reproductive  rate  depends  on  the  time  index  in  addition  to  the  NAO  index,  and 
the  slope  parameter  associated  with  the  time  index  in  these  models  is  negative.  Therefore, 
the  data  also  suggest  that  the  reproductive  rate  has  been  declining  over  time  independent 
of  the  effect  of  NAO.  Finally,  these  models  have  survival  probability  of  mothers  that  depend 
on  NAO;  there  is  evidence  that  environmental  forcing  is  acting  on  the  survival  probability 
of  mothers. 

The  AAIC  values  for  the  first  three  models  in  Table  6.3.1  are  less  than  2;  so  there  is 
not  strong  evidence  to  exclude  these  models  from  potentially  the  best  model  in  the  data 
(Burnham  &  Anderson,  1998).  Therefore,  the  rest  of  our  analyses  focuses  on  these  models 
separately  and  compares  the  results  to  find  the  influence  of  model  selection  uncertainty.  We 
denote  these  three  statistical  models  by  M\,  M2,  and  M3  as  indicated  in  Table  6.3.1. 

6.3.2  Estimated  probabilities 
Capture  probability  model 

Figure  6.3  shows  the  estimated  stage-specific  capture  probabilities  of  the  three  best  models; 
they  are  almost  identical.  This  suggests  that  the  capture  probability  estimates  are  not 
affected  by  the  differences  in  transition  probability  models.  The  capture  probability  of  all 
stages  appears  to  have  increased  over  the  course  of  the  study.  This  reflects  the  fact  that 
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the  sampling  effort  has  increased  (see  Chapter  5).  Compared  with  other  stages,  the  capture 
probability  of  mothers  is  very  high  and  close  to  1  (except  in  some  early  years).  This  indicates 
that  mothers  are  easier  to  observe  than  other  stages,  or  the  sampling  program  is  particularly 
effective  for  detecting  them.  The  capture  probabilities  of  mature  females  during  the  year 
immediately  after  the  mother  stage  and  other  years  (stage  5  and  3,  respectively)  appear  to 
be  different.  A  likelihood  ratio  test  indicates  that  this  difference  is  significant.  Individuals 
in  stage  3  are  available  for  reproduction  in  the  following  year,  but  those  in  stage  5  are  not 
because  of  reproduction  in  the  previous  year.  This  difference  could  influence  the  behavior 
of  females  and  affect  the  capturability  of  the  individuals.  All  mature  females  (stage  3  and 
5)  are  assumed  to  have  the  same  capture  probability  in  the  previous  study  in  Chapter  5; 
this  difference  in  the  capture  probability  could  have  caused  slight  underestimation  of  the 
previous  estimate  of  survival  probability  of  mothers  (stage  4). 

Transition  probability  models 

Figure  6.4  shows  estimated  stage-specific  survival  probabilities.  The  survival  probabilities  of 
calves  are  different  among  the  three  models.  In  the  first  model  (Mi),  the  survival  probability 
is  negatively  correlated  with  the  NAO  index.  Because  the  index  has  a  slightly  increasing 
trend  over  time,  the  survival  probability  also  has  a  slightly  declining  trend.  In  the  second 
model  (M2),  the  survival  probability  has  a  monotonically  declining  trend.  In  the  third 
model  (M3),  the  calf  survival  is  a  function  of  NAO  only.  We  do  not  have  enough  data 
to  test  whether  the  survival  probability  of  calves  is  significantly  correlated  with  the  NAO 
index  or  not;  however,  the  three  models  suggest  that  there  is  at  least  a  slightly  declining 
trend. 

The  survival  probabilities  of  mothers  in  the  best  (Mi)  and  the  second  best  (M2)  models 
decline  monotonically  with  time.  The  same  trend  was  also  found  in  Chapter  5.  On  the  other 
hand,  the  survival  probability  of  mothers  in  the  third  model  (M3)  is  negatively  correlated 
with  both  the  time  and  NAO  indices.  This  model  shows  that  the  survival  probability  in 
1996  was  dramatically  low.  The  fact  that  there  were  unusually  high  observed  mortalities  in 
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Figure  6.4:  Stage-specific  survival  probabilities:  (a)  female  under  Mi,  (b)  female  under  M2, 
(c)  female  under  M3,  (d)  male  under  Mi,  (e)  male  under  M2,  and  (f)  male  under  M3.  The 
curves  are  for  calves  (squares),  juvenile  (triangles),  mature  (circles),  mother  (diamonds), 
and  year  after  the  mother  stage  (stars). 
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the  winter  of  1996  (Waring  et  al.,  1999)  indicates  that  model  M3  may  be  more  accurately 
reflecting  the  true  state  of  the  population.  In  all  three  models,  the  survival  probability  of 
mothers  is  declining.  This  decline  is  significant  because  the  slope  parameters  associated 
with  time  (72  in  equation  (6.10))  in  all  three  models  are  significantly  less  than  0  based  on 
the  likelihood  ratio  tests. 

Figure  6.5  shows  estimated  conditional  transition  probabilities  from  stage  3  to  stage 
4  of  the  three  best  models.  The  conditional  transition  probability  gives  the  probability  of 
mature  females  becoming  mother  provided  they  survive  until  the  following  year.  In  all  three 
models,  the  probability  is  correlated  with  both  the  NAO  and  time  indices.  This  result  is 
consistent  with  the  previous  studies  by  Knowlton  et  al.  (1994);  Pershing  &  Greene  (2002). 
This  suggests  there  is  declining  trend  in  the  reproductive  rate  upon  which  the  effect  of 
climate  fluctuations  are  superimposed. 

6.3.3  Uncertainty  in  parameter  estimates 

The  estimated  probabilities  in  the  previous  section  have  associated  uncertainties  that  we 
would  like  evaluate  in  terms  of  some  indicator  of  population  health.  With  a  lack  of  an 
universal  indicator,  we  use  the  asymptotic  annual  population  growth  rate  (A)  as  a  proxy 
for  population  health. 

We  drew  the  parameters  associated  with  the  stage-specific  survival  probabilities  and 
conditional  transition  probabilities  from  a  multivariate  normal  distribution.  The  covari¬ 
ance  matrix  for  the  distribution  was  approximated  by  the  inverse  of  the  Hessian  matrix 
calculated  at  the  last  iteration  in  the  numerical  maximum  likelihood  method  (see  Burnham 
et  al.,  1987).  The  simulated  parameters  were  transformed  into  the  stage-specific  survival 
probabilities  and  the  conditional  transition  probabilities  setting  the  NAO  index  to  the  mean 
value  between  1980  to  1997  and  the  time  to  be  1  (i.e.  the  level  in  1997)  in  equations  (6.9)- 
(6.11).  By  setting  the  NAO  index  to  the  mean  value,  we  removed  the  variability  associated 
with  the  index.  By  setting  the  time  to  the  1997  level,  we  assume  that  mean  environmental 
conditions  remain  at  1997  levels  in  future  years.  With  these  probabilities,  we  constructed 
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a  population  projection  matrix.  From  the  matrix,  we  calculated  an  asymptotic  population 
growth  rate  (A),  which  is  given  by  the  dominant  eigenvalue  of  the  projection  matrix  (see 
Caswell,  2001).  This  was  repeated  1000  times  for  each  of  the  three  models  selected  in  the 
previous  section. 

The  result  shows  that  there  are  79%  (Mi),  86%  (M2),  and  81%  (M3)  chances  that  the 
asymptotic  population  growth  rate  (A)  is  less  than  1 .  When  A  was  less  than  1 ,  we  calculated 
the  quasi-extinction  time  (r<f).  Here,  we  define  the  quasi-extinction  time  to  be  the  time  that 
it  takes  a  population  to  be  reduced  from  150  females  distributed  according  to  stable  stage 
distribution  to  a  single  female.  The  quasi-extinction  time  (r^)  was  calculated  as 

fd  =  (6.12) 

log  A 

Figure  6.6  shows  the  distribution  of  the  quasi-extinction  time  (?d)  conditional  on  A  <  1. 
This  figure  shows  that  there  is  a  wide  range  in  fj,.  It  also  shows  that  when  the  population 
goes  extinct,  it  disappears  relatively  quickly.  In  fact,  of  the  1000  simulations  (including  the 
cases  in  which  A  >  1  and  A  <  1),  there  is  50%  chance  of  the  population  going  extinct  within 
236  (Mi),  209  (M2),  and  176  (Mb)  years. 

6.3.4  Environmental  Stochasticity 

The  quasi-extinction  time  calculation  in  the  previous  section  assumes  a  deterministic  change 
in  the  population  size  and  does  not  include  any  environmental  stochasticity.  However, 
populations  live  in  stochastic  environments,  and  we  would  like  to  assess  the  effect  of  this 
stochasticity  on  population  health.  Here,  we  analyze  the  effect  of  adding  variability  to  the 
mean  stage-specific  survival  and  conditional  transition  probabilities  in  the  previous  section. 

The  variability  in  the  probabilities  (both  stage-specific  survival  probabilities  and  the 
conditional  transition  probabilities)  due  to  environmental  stochasticity  was  approximated 
by  the  deviations  in  saturated  models,  in  which  probabilities  change  freely  from  one  year 
to  the  next,  from  one  of  the  three  best  models  (Mi,  M2,  and  M3).  First,  we  calculated  the 
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Time  to  Quasi-extinction  (Year) 


Figure  6.6:  Distributions  of  quasi-extinction  time  showing  parameter  uncertainty  under 
Mi  (a),  M2  (b),  and  M3  (c).  Parameters  generated  with  estimated  means  and  associated 
covariance  matrix  assuming  their  multivariate  Normal  distribution. 
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probabilities  using  the  parameters  estimated  with  the  saturated  model.  We  also  calculated 
the  probabilities  using  the  parameters  estimated  with  one  of  the  three  models  using  the 
mean  NAO  index  but  the  time  index  changing  with  time  in  equations  (6.9)-(6.11).  Then  we 
subtracted  these  probabilities  from  the  corresponding  probabilities  of  the  saturated  model  to 
calculate  the  deviations.  Finally,  these  deviations  were  added  to  the  probabilities  calculated 
using  parameters  of  the  same  model  (i.e.  one  of  the  three  best  models)  but  setting  the  NAO 
and  time  indices  to  the  mean  and  the  1997  level,  respectively.  This  resulted  in  17  values  for 
each  probability.  Two  of  the  values  associated  with  £34  in  each  model  were  negative;  they 
were  set  to  0.  Using  these  probabilities,  we  constructed  17  population  projection  matrices. 

The  17  matrices  approximate  the  variability  in  a  population  projection  matrix  around 
the  matrix  that  would  have  been  expected  in  the  environmental  condition  of  1997  but 
under  the  average  climate  level.  It  should  be  noted  that  we  did  not  subtract  the  sampling 
variations  from  the  estimates.  Therefore,  the  17  matrices  are  over-estimating  the  variability 
due  to  environmental  stochasticity. 

Using  the  17  population  projection  matrices,  we  calculated  the  cumulative  probability 
density  for  time  to  quasi-extinction  assuming  the  changes  in  log  of  the  total  population  size 
are  approximated  by  a  diffusion  process  (e.g.  Caswell,  2001;  Dennis  et  al.,  1991).  For  this 
calculation,  the  stochastic  growth  rate  was  calculated  from  a  population  growth  rate  of  the 
average  matrix  using  the  approximation  developed  by  Tuljapurkar  (1990)  and  the  diffusion 
coefficient  (a2)  was  approximated  by 

a2  =  2(logp-  log  As),  (6.13) 

where  log  p  is  the  population  growth  rate  of  the  average  matrix  and  log  As  is  the  stochastic 
population  growth  rate. 

The  mean  deviations  added  to  all  stage-specific  survival  probabilities  in  all  three  models 
except  the  one  for  survival  probability  of  mothers  in  M3  were  negative.  Therefore,  by 
adding  the  deviations,  we  reduced  the  survival  probabilities  on  average.  Paradoxically, 
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the  deterministic  growth  rate  of  the  projection  matrix  averaged  over  the  17  matrices  was 
higher  than  the  deterministic  growth  rate  of  the  matrix  before  adding  the  environmental 
variability.  As  a  result,  the  quasi-extinction  time  calculated  using  equation  (6.12)  increased 
from  361  based  on  the  mean  estimates  for  1997  to  447  based  on  the  average  of  the  17 
matrices  constructed  using  Mi  as  the  trend  model,  from  354  based  on  the  mean  estimates 
for  1997  to  437  based  on  the  average  of  the  17  matrices  constructed  using  M2  as  the  trend 
model,  and  from  196  based  on  the  mean  estimates  for  1997  to  248  based  on  the  average 
of  the  17  matrices  constructed  using  M3  as  the  trend  model.  This  effect  comes  from  the 
skewed  distributions  of  the  survival  probabilities  and  the  negative  mean  deviations  added 
to  <43.  When  443  is  reduced,  fewer  individuals  make  the  transition  into  the  mother  stage  in 
which  individuals  face  higher  risk  of  mortality. 

Figure  6.7  shows  the  cumulative  probability  density  for  quasi-extinction.  In  each  figure, 
the  vertical  line  on  the  left  indicates  the  median  time  to  the  quasi-extinction  and  the  line 
on  the  right  indicates  the  quasi-extinction  time  calculated  using  the  deterministic  model 
(equation  (6.12)).  These  two  values  are  close  to  each  other  in  all  models. 

6.3.5  Demographic  stochasticity 

The  extinction  calculations  so  far  do  not  include  the  effect  of  the  small  population  size 
(i.e.  demographic  stochasticity).  In  this  section,  we  analyze  this  effect  on  the  population 
viability.  Changes  in  population  size  are  determined  by  the  fates  of  individuals;  those  fates 
are  determined  stochastically.  For  example,  with  some  probability  each  individual  survives 
or  die,  and  with  some  probability  an  individual  makes  a  transition  into  a  different  stage  or 
remains  in  the  same  stage;  these  stochasticities  occur  even  if  the  environment  is  constant. 
When  population  size  is  large,  the  effect  of  these  events  are  averaged  over  many  individuals, 
and  as  a  result,  the  population  size  changes  smoothly.  However,  when  the  population  size 
is  small,  the  fate  of  each  individual  has  a  magnified  effect  on  the  total  population  size.  This 
is  called  demographic  stochasticity. 

The  effect  of  demographic  stochasticity  on  the  viability  of  the  North  Atlantic  right 
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Figure  6.7:  Cumulative  probability  density  for  quasi-extinction  time  showing  the  effect 
of  environmental  stochasticity  under  M\  (a),  M2  (b),  and  M3  (c).  The  left  vertical  line 
shows  median,  and  the  right  vertical  line  shows  quasi-extinction  time  calculated  under  the 
deterministic  model,  (see  text  for  the  method) 
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whales  has  been  analyzed  in  Chapter  5.  Here,  we  repeat  the  same  analyses  based  on  the 
new  estimates.  For  this  analysis,  we  used  the  average  matrix  used  in  the  previous  section. 
These  estimates  were  analyzed  using  the  multi- type  branching  process  (see  Caswell,  2001). 

Figure  6.8  shows  the  cumulative  probability  density  for  the  extinction  time.  As  in  the 
previous  section,  the  vertical  line  on  the  left  is  the  median  time  to  extinction,  and  the 
vertical  line  on  the  right  is  the  quasi-extinction  time  based  on  the  deterministic  calculation 
(equation  6.12).  It  should  be  noted  that  we  are  calculating  the  time  to  extinction  (not  quasi¬ 
extinction)  using  the  branching  process  calculation.  The  extinction  time  is  defined  here  as 
the  time  it  takes  a  population  consisting  of  150  females  distributed  according  to  stable  stage 
distribution  to  completely  disappear.  Even  the  comparison  of  the  quasi-extinction  and  the 
extinction  times,  it  is  clear  that  there  is  a  larger  effect  of  the  demographic  stochasticity  on 
the  population  viability  than  the  environmental  stochasticity. 

6.3.6  Environmental  and  demographic  stochasticity 

After  the  separate  analyses  of  environmental  and  demographic  stochasticities,  the  next 
obvious  step  is  to  analyze  the  combined  effects  on  the  population  viability.  When  both 
environmental  stochasticities  and  demographic  stochasticities  are  combined,  there  is  no 
practical  analytical  method  to  calculate  the  extinction  time.  So  we  rely  on  simulations. 

We  started  the  population  with  150  females  distributed  according  to  the  stable  stage 
distribution  of  the  average  matrix  in  the  previous  section.  Then  for  each  year,  we  randomly 
selected  one  of  the  17  matrices  with  equal  probabilities.  Using  the  selected  matrix,  we 
projected  the  number  of  individuals  in  the  zth  stage  in  time  t  +  1  coming  from  stage  j 
in  year  t  assuming  multinomial  distributions  with  the  number  of  individuals  in  stage  j 
in  year  t  and  probabilities  faj  as  parameters.  This  was  repeated  for  all  stages  for  each 
year  of  projection.  The  number  of  calves  in  year  t  +  1  was  calculated  using  the  binomial 
distribution  with  the  number  of  individual  in  mature  stage  (stage  3)  in  year  y  and  F3  as 
parameters.  The  population  projection  was  repeated  until  the  population  goes  to  quasi- 
and  actual  extinction,  and  the  years  of  these  extinction  events  were  noted.  This  process 
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Figure  6.8:  Cumulative  probability  density  for  quasi-extinction  time  showing  the  effect 
of  demographic  stochasticity  under  M\  (a),  M2  (b),  and  M3  (c).  The  left  vertical  line 
shows  median,  and  the  right  vertical  line  shows  quasi-extinction  time  calculated  under  the 
deterministic  model,  (see  text  for  the  method) 
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was  repeated  for  500  times  using  the  parameter  estimates  of  each  model. 

Figure  6.9  shows  50  randomly  selected  simulation  paths  of  the  total  population  size. 
The  solid  lines  in  the  figures  indicate  the  deterministic  path  of  the  population  size.  Figure 
6.10  shows  the  distributions  of  actual  extinction  time.  As  expected,  the  combination  of 
demographic  and  environmental  stochasticities  creates  a  range  of  projected  extinction  times. 

The  above  simulation  assumes  that  mature  females  giving  birth  and  making  transi¬ 
tion  into  the  mother  stage  are  independent.  This  assumption  was  also  made  in  all  of  the 
calculations  so  far.  However,  in  reality,  every  individual  that  makes  the  transition  from 
mature  to  mother  stage  also  gives  birth.  Thus  the  reproduction  and  transition  are  corre¬ 
lated.  Although  incorporating  this  correlation  in  the  diffusion  approximation  model  and 
the  branching  process  calculation  is  difficult,  it  can  be  incorporated  easily  in  a  simulation. 
Therefore,  we  repeat  the  same  simulation  as  before  but  calculating  the  number  of  calves 
in  year  y  4-  1  differently  to  account  for  the  correlation.  In  this  calculation,  when  individ¬ 
uals  made  the  transition  from  mature  (stage  3)  to  mother  (stage  4)  stages,  we  used  the 
binomial  distribution  with  number  of  individuals  that  made  the  transition  and  0.5</>54°  5  as 
parameters  to  decide  how  many  individuals  appear  as  calves  in  the  following  year. 

Assuming  the  independence  between  the  reproduction  and  the  transition,  the  median 
times  to  quasi-extinction  were  317  {Mi),  315  (M2),  and  194  (M3)  years,  and  the  median 
times  to  actual  extinction  were  355  (Mi),  352  (M2),  and  217  (M3)  years.  However,  when 
we  assume  the  correlation,  the  median  times  to  quasi-extinction  increase  to  338  (Mi),  333 
(M2),  and  201  (M3)  years,  and  the  median  time  to  the  extinction  increase  to  380  (Mi), 
377  (M2),  and  229  (M3)  years.  These  increases  result  from  the  fact  that  individuals  in 
the  mother  stage  have  very  low  survival  probability  compared  with  other  stages.  If  the 
probability  of  individuals  risking  into  the  mother  stage  and  reproduction  are  correlated, 
there  is  a  less  chance  that  mothers  die  without  giving  birth  and  a  higher  chance  that  the 
population  size  does  not  decline.  We  have  not  included  the  correlation  between  survival  of 
mothers  and  calves  in  our  model.  The  right  whale  mothers  nurse  their  calves  for  almost  a 
year.  It  is  very  likely  that  calves  die  if  mothers  die.  This  correlation  would  have  an  opposite 
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Figure  6.9:  Stochastic  simulations  of  the  total  population  size  with  both  environmental  and 
demographic  stochasticity  under  Mi  (a),  M2  (b),  and  M3  (c).  50  randomly  selected  paths 
are  shown.  The  simulation  does  not  assume  the  correlation  between  transition  into  the 
mother  stage  and  fertility.  The  solid  line  is  a  path  of  the  deterministic  model. 
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Figure  6.10:  Distribution  of  quasi-extinction  time  showing  the  effect  of  both  environmental 
and  demographic  stochasticity  under  M\  (a),  M2  (b),  and  M3  (c).  The  simulation  does  not 
assume  the  correlation  between  transition  into  the  mother  stage  and  fertility. 
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effect  as  the  other  correlation  and  reduce  the  time  to  extinction. 


6.3.7  More  simulations 

We  used  three  different  analytical  methods  to  calculate  expected  time  to  extinction  (i.e. 
deterministic,  diffusion  approximation,  and  branching  process);  however,  the  comparisons 
of  the  results  are  not  trivial.  For  example,  the  deterministic  and  diffusion  approximation 
methods  calculate  the  time  to  quasi-extinction,  but  the  branching  process  method  calculates 
the  time  to  actual  extinction.  Furthermore,  the  three  methods  may  have  different  implicit 
assumptions  or  use  approximations  that  may  be  sensitive  to  different  factors.  Theoretical 
studies  of  these  methods  are  beyond  a  scope  of  this  study.  Nevertheless,  it  would  be  nice 
if  we  could  compare  the  results  based  on  the  same  method.  Fortunately,  the  simulation 
method  used  in  the  previous  section  also  allows  us  to  analyze  environmental  or  demographic 
stochasticity  alone. 

For  the  analysis  of  the  demographic  stochasticity  alone,  we  used  the  method  in  the 
previous  section  but  only  using  the  average  matrix  alone.  Here  we  assumed  no  correlation 
between  the  reproduction  and  the  transition  from  mature  to  mother  stage  to  be  consistent 
with  the  branching  process  calculation.  Using  this  method,  we  calculated  actual  and  quasi¬ 
extinction  time. 

For  the  analysis  of  the  environmental  stochasticity  alone,  we  randomly  selected  one  of  the 
17  matrices  in  the  previous  section  for  each  year.  However,  instead  of  using  the  multinomial 
distribution  to  project  the  population  in  the  following  year,  we  multiplied  the  projection 
matrix  to  a  vector  containing  number  of  individuals  in  each  stage  (equation  (6.2)).  For 
the  initial  population,  we  used  150  individuals  distributed  according  to  the  stable  stage 
distribution  of  the  average  matrix.  We  repeated  this  process  until  the  sum  of  the  total 
population  size  becomes  1.  We  noted  this  year  as  a  quasi-extinction  time.  This  process  was 
repeated  500  times  for  each  model. 

Table  6.3.7  summarizes  all  the  extinction  time  calculations  made  in  this  paper.  The 
result  shows  that  the  analytical  calculations  using  diffusion  approximation  and  branching 
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Table  6.2:  Results  of  extinction  time  calculations 


Method 

Type  of 

Factors 

Median  Years  to  Extinction 

Extinction 

Analyzed 

Mi 

m2 

m3 

Deterministic 

Quasi 

447 

438 

248 

Diffusion  Process 

Quasi 

ES 

438 

429 

245 

Branching  Process 

Actual 

DS 

359 

355 

219 

Simulation 

Quasi 

ES 

437 

430 

250 

Simulation 

Actual 

DS 

359 

357 

224 

Simulation 

Quasi 

DS 

313 

318 

199 

Simulation 

Quasi 

DS+ES 

338 

333 

201 

Simulation 

Actual 

DS+ES 

380 

377 

229 

Simulation 

Quasi 

DS+ES+CR 

317 

315 

194 

Simulation 

Actual 

DS+ES+CR 

355 

352 

217 

DS:  demographic  stochasticity. 

ED:  environmental  stochasticity. 

CR:  correlation  between  reproduction  and  transition  into  mother  stage. 


process  methods  agree  very  well  with  the  simulations.  The  comparisons  of  actual  and  quasi¬ 
extinction  times  for  the  demographic  stochasticity  alone  indicate  that  there  is  about  20  to 
40  years  difference  between  them. 

6.4  Discussion 

We  analyzed  the  North  Atlantic  right  whale  data  using  64  models  that  differ  in  hypotheses 
on  how  reproductive  rate  and  survival  probability  of  calves  and  mothers  changed  over  time. 
Based  on  AIC,  we  selected  the  three  statistical  models  that  were  supported  by  the  data 
best.  According  to  these  models,  the  survival  probability  of  mothers  has  declined  over 
time  and  the  reproductive  rate  of  mature  females  has  been  correlated  with  the  climate 
index.  The  three  models  disagree  on  whether  the  survival  probability  of  mothers  has  been 
correlated  with  the  climate  index  and  also  disagree  on  how  the  survival  probability  of  calves 
has  changed  over  time.  However,  throughout  the  series  of  the  population  viability  analyses 
in  this  paper,  results  based  on  the  three  models  did  not  give  qualitatively  different  answers. 
Based  on  the  three  models,  we  analyzed  the  uncertainties  in  parameter  estimates  using  a 
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parametric  bootstrap  method.  For  this  purpose,  we  used  an  asymptotic  population  growth 
rate  (A)  as  an  index.  According  to  the  analysis,  all  three  models  agree  on  a  80%  chance 
of  extinction  (A  <  1)  provided  that  climatic  conditions  are  maintained  at  the  mean  level 
between  1980  and  1997  and  that  other  environmental  conditions  are  maintained  at  the  1997 
level  in  future  years. 

Using  the  estimates  based  on  the  three  models,  we  analyzed  the  effects  of  environmental 
and  demographic  stochasticity  on  the  viability  of  the  population.  We  used  four  different 
methods  to  calculate  the  quasi-extinction  and  actual  extinction  times  (Table  6.3.7).  The 
results  reveal  that  there  is  little  effect  of  environmental  stochasticity  on  the  median  time 
to  quasi-extinction.  The  environmental  stochasticity  appears  only  to  increase  the  range  of 
quasi-extinction  time  projections.  Considering  that  North  Atlantic  right  whales  are  long- 
lived  animals,  and  we  are  only  adding  noise  at  a  small  time-scale,  this  result  is  reasonable. 
On  the  other  hand,  demographic  stochasticity  has  much  larger  effects  on  population  vi¬ 
ability.  Simulation  results  suggest  that  there  is  an  80  to  90%  chance  that  demographic 
stochasticity  reduces  the  time  to  quasi-extinction,  and  there  is  a  20  to  30%  reduction  in 
median  time  to  quasi-extinction. 

The  North  Atlantic  right  whale  is  seriously  endangered  not  only  because  of  its  small 
population  size,  but  also  because  of  the  trend  in  the  population  growth  rate  (Caswell  et  al., 
1999,  Chapter  5).  Previous  studies  suggested  that  the  population  is  now  slowly  declining. 
In  the  current  study,  we  analyzed  the  data  on  this  population  taking  additional  factors 
into  consideration.  Those  additional  factors  are  model  selection  uncertainty,  the  effect  of 
environmental  stochasticity,  and  some  heterogeneity  in  the  capture  probability  of  mature 
females.  With  all  these  factors  considered,  the  results  do  not  suggest  a  bright  future  for 
this  population. 

Recently  a  few  major  management  actions  have  been  implemented.  They  include  a 
disentanglement  network  system,  a  mandatory  ship  reporting  system,  and  a  temporary 
closure  of  areas  from  fishing  activities  (Waring  et  al.,  1999).  These  management  actions 
were  implemented  in  the  mid  to  late  1990's,  and  their  effects  will  not  be  detectable  for  a 
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few  more  years.  Saving  only  a  few  whales  will  protect  the  population  (Chapter  5).  We  hope 
that  these  management  actions  are  effective  and  that,  in  the  future,  we  will  be  reporting 
the  effects  of  successful  management  on  the  viability  of  the  North  Atlantic  right  whales. 


Chapter  7 

Conclusion 


The  main  chapters  (Chapters  2-6)  were  written  for  possible  publications  in  journals,  and 
they  contain  conclusions  of  their  own.  Therefore,  I  present  a  brief  overview  of  this  thesis  in 
this  chapter. 

In  this  thesis,  I  presented  new  methods  to  analyze  mark-recapture  data  (Chapters  2,  3, 
and  4)  and  applications  of  the  methods  to  actual  data  (Chapters  5  and  6).  These  methods 
link  mark-recapture  data  and  existing  population  models  that  allow  modern  demographic 
analyses  (e.g.  Caswell,  2001).  I  presented  selected  demographic  analyses  such  as  sensitiv¬ 
ity,  branching  process,  and  environmental  stochasticity  analyses  based  on  the  estimated 
parameters  of  the  right  whale  models  in  Chapter  5  and  6.  The  study  of  the  right  whale 
population  presented  in  Chapter  5  revealed  some  interesting  facts  about  the  status  of  the 
population.  One  of  the  most  interesting  findings  is  that  the  declining  survival  probability 
found  by  Caswell  et  al.  (1999)  appears  to  be  caused  by  the  increased  mortality  rate  of 
females  that  are  nursing  their  offspring.  I  believe  that  this  and  other  results  in  the  two 
chapters  demonstrate  the  usefulness  of  the  mark-recapture  methods  presented  in  the  earlier 
chapters. 

As  marking  techniques  improve,  individual  mark-recapture  data  are  becoming  very  com¬ 
mon.  Such  data  are  already  available  for  wandering  albatross  ( Diomedia  exulans )  and  Am¬ 
sterdam  albatross  ( D .  amsterdamensis )  (Weimerskirch  et  al.,  1997),  Atlantic  humpback 
whales  (Buckland,  1990),  kittiwake  gull  ( Rissa  tridactyla )  (Aebischer  &  Coulson,  1990), 
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Canada  geese  ( Branta  canadensis )  (Hestbeck  et  al.,  1991),  tundra  swan  ( Cygnus  columbi- 
nus )  (Nichols  et  al.,  1992a),  roseate  tern  ( Sterna  dougallii )  (Spendelow  et  al.,  1995),  black¬ 
headed  gull  ( Larus  ridibundus )  (Prevot-Julliard  et  al.,  1998),  Florida  manatee  ( Trichechus 
manatus  latirostris)  (Langtimm  et  al.,  1998),  pig  frog  ( Rana  grylio)  (Wood  et  al.,  1998), 
Weddell  seal  ( Leptonychotes  weddellii )  (Hastings  &  Testa,  1998),  grizzly  bear  ( Ursus  arc- 
tos)  (Pease  &  Mattson,  1999),  grey  seal  ( Halichoerus  grypus )  (Schwarz  &  Stobo,  2000),  and 
northern  spotted  owl  ( Strix  occidentalis  caurina )  (Franklin  et  al.,  2000)  to  name  some  exam¬ 
ples.  The  methods  I  presented  allow  reanalysis  of  these  existing  data  taking  full  advantage 
of  the  demographic  analyses  and  also  allow  planning  of  future  sampling. 

I  here  show  that  mark-recapture  data  contain  a  great  deal  of  information.  They  are 
one  of  the  best  data  to  collect  when  studying  a  population.  I  was  able  to  find  a  change 
in  survival  probability  in  only  one  fife  stage  of  individuals  in  a  population  that  consists  of 
about  300  individuals.  This  discovery  would  not  have  been  possible  if  there  were  not  the 
mark-recapture  data  and  the  appropriate  analysis  method  presented  in  this  thesis.  However, 
mark-recapture  data  are  not  always  available.  One  obvious  limitation  of  the  method  is  that 
individuals  must  be  uniquely  marked.  When  individuals  cannot  be  marked,  an  alternative 
type  of  data  is  count  data,  which  record  how  many  individuals  in  each  stage  were  observed 
at  each  sampling  occasion.  Recently,  Zeng  et  al.  (1998)  and  DeValpine  &  Hastings  (2002) 
developed  a  method  to  estimate  population  parameters  by  fitting  state-space  models  to 
count  data.  Because  count  data  are  easier  to  collect  and  more  available  than  mark-recapture 
data,  I  suspect  their  method  will  become  very  important  in  ecological  studies.  However,  I 
am  not  sure  which  type  of  data  are  superior  given  the  same  amount  of  sampling  effort.  I 
think  the  answer  will  depend  on  the  size  of  a  population  and  detectability  of  individuals  (i.e. 
capture  probability).  It  is  also  conceivable  to  collect  both  mark-recapture  and  count  data 
on  the  same  population.  Then  a  single  likelihood  function  for  both  data  can  be  expressed 
to  estimate  parameters.  Under  such  circumstance,  I  am  interested  to  find  what  is  the  best 
allocation  of  effort  between  the  two  types  of  data. 

Collecting  data  in  a  population  study,  whether  they  are  mark-recapture  data  or  count 
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data,  is  often  an  expensive  process.  The  cost  of  the  collecting  data  should  warrant  the  most 
sophisticated  analyses  available.  Once  data  are  collected,  they  should  be  analyzed  to  a  full 
extent.  Unfortunately,  ecological  data  are  often  under-utilized.  In  my  opinion,  this  comes 
from  the  lack  of  emphasis  in  studies  that  try  to  link  data  and  models.  This  is  the  area  in 
which  I  am  hoping  to  make  contributions  as  an  ecologist/biological  oceanographer.  This 
thesis  is  the  first  of  such  contributions. 
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Mark-recapture  analyses  of  populations  are  an  important  tool  in  population  biology.  In  this  thesis,  I  extend  mark-recapture 
analyses  to  provide  a  link  between  mark-recapture  data  and  demographic  models  such  as  matrix  population  models  and 
integrodifference  equation  models.  I  resolve  some  issues  that  are  commonly  encountered  during  sampling,  such  as  the  fact 
that  the  sex  or  life-stage  of  some  individuals  is  unknown  during  some  sampling  occasions  and  that  individuals  become 
unobservable  during  some  life  stages.  I  introduce  a  stage  structure  that  permits  simple  conversion  of  estimated  transition 
probabilities  into  a  matrix  population  model.  I  describe  a  simple  algorithm  to  simplify  programming  for  parameter 
estimation.  I  also  introduce  a  method  to  estimate  the  distribution  of  dispersal  displacements  (dispersal  kernel)  from 
mark-recapture  data.  I  use  some  of  the  above  methods  to  estimate  the  vital  rates  of  the  North  Atlantic  right  whale 
(. Eubalaena  glacialis).  To  the  estimated  vital  rates,  I  apply  demographic  analyses  including  population  viability  analyses 
and  sensitivity  analysis.  Finally,  I  compare  effects  of  environmental  and  demographic  stochasticities  on  the  viability  of 
the  population. 


17.  Document  Analysis  a.  Descriptors 
Mark-recapture  statistics 
demographic  analysis 
matrix  population  models 

b.  Identifiers/Open-Ended  Terms 


c.  COSATI  Field/Group 


18.  Availability  Statement 

Approved  for  publication;  distribution  unlimited. 


(See  ANSI-Z39.18)  £ 


19.  Security  Class  (This  Report) 

UNCLASSIFIED 

21.  No.  of  Pages 

138 

20.  Security  Class  (This  Page) 

22.  Price 

See  Instructions  on  Reverse 


OPTIONAL  FORM  272  (4-77) 

(Formerly  NTIS-35) 
Department  of  Commerce 


